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Abstract 

The paper presents a detailed review of the smooth particle hydrodynamics (SPH) 
method with particular focus on its astrophysical applications. We start by intro- 
ducing the basic ideas and concepts and thereby outline all ingredients that are nec- 
essary for a practical implementation of the method in a working SPH code. Much 
of SPH's success relies on its excellent conservation properties and therefore the nu- 
merical conservation of physical invariants receives much attention throughout this 
review. The self-consistent derivation of the SPH equations from the Lagrangian 
of an ideal fluid is the common theme of the remainder of the text. We derive a 
modern, Newtonian SPH formulation from the Lagrangian of an ideal fluid. It ac- 
counts for changes of the local resolution lengths which result in corrective, so-called 
"grad-h-terms". We extend this strategy to special relativity for which we derive 
the corresponding grad-h equation set. The variational approach is further applied 
to a general-relativistic fluid evolving in a fixed, curved background space-time. 
Particular care is taken to explicitely derive all relevant equations in a coherent 
way. 
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1 Introduction 



Much of what we observe in the physical Universe has been shaped by fluid 
dynamical processes. Prom the hot gas in galaxy clusters to the internal struc- 
tures of their consitucnt galaxies down to their stars and planets, all has been 
formed by the interplay between gravity, gas dynamics and further physical 
processes such as the interaction with radiation, nuclear burning or magnetic 
flelds. The latter processes often involve intrinsic length and time scales that 
are dramatically different from those of the gas dynamical processes, therefore 
many astrophysical problems are prime examples of multi-scale and multi- 
physics challenges. The complexity of the involved physical processes and the 
lack of symmetry usually prohibit analytical treatments and only numerical 
approaches are feasible. 

Although fluid dynamics is also crucial for many technical applications, their 

requirements usually differ substantially from those of astrophysics and this 
also enters the design of the numerical methods. Typical astrophysical require- 
ments include: 

• Since flxed boundaries are usually absent, flow geometries are determined 
by the interplay between different physical processes such as gas dynamics 
and (self-)gravity which often lead to complicated, dynamically changing 
flow geometries. Thus, a high spatial adaptivity is often required from as- 
trophysical hydodynamics schemes. 

• Shocks often crucially determine the evolution of cosmic objects. Examples 
include in supernova remnants or the Earth's magnetosphere. 

• Physical quantities can vary by many orders of magnitudes between dif- 
ferent regions of the simulation domain. This requires a particularly high 
robustness of the numerical scheme. 

• In many astrophysical problems the numerical conservation of physically 
conserved quantities determines the success and the reliability of a computer 
simulation. Consider, for example, a molecular gas cloud that collapses un- 
der the influence of its own gravity to form stars. If the simulation for some 
reason dissipates angular momentum, a collapsing, self-gravitating portion 
of gas may form just a single stellar object instead of a multiple system of 
stars and it will thus produce a qualitatively wrong result. 

• Many astrophysical questions require dealing with physical processes beyond 
gas dynamics and self-gravity. A physically intuitive and flexible formula- 
tion of the numerics can substantially facilitate the implementation of new 
physics modules into existing codes. 

No numerical method performs equally well at each of the above requirements, 
therefore, the choice of the best-suited numerical approach can often save a 
tremendous amount of effort in obtaining reliable results. Therefore: horses 
for courses. 
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In the following, the smooth particle hydrodynamics (SPH) method [1-6], a 
completely mesh-free approach to solve the hydrodynamic equations is dis- 
cussed in detail. Its conservation properties are a major strength of SPH, 
therefore "hard-wired" conservation receives much attention throughout this 
text. The derivation of the SPH equations (in the absence of dissipation) re- 
quires nothing more than a suitable Lagrangian, a density prescription that 
depends on the coordinates and the first law of thermodynamics. The resulting 
equations conserve the physically conserved quantities even in their discretized 
form, provided that the original Lagrangian possessed the correct symmetries. 
Therefore, derivations from Lagrangians play a central role in our discussion 
of the subject. 

This review has emerged from a lecture series on "Computational relativistic 
astrophysics" as part of a Doctoral Training Programme on the "Physics of 
Compact Stars" that was held in summer 2007 at the European Centre for 
Theoretical Studies in Nuclear Physics and Related Areas (ECT*) in Trento, 
Italy. In the pedagogical spirit of this lecture scries the text is kept in "lec- 
ture" rather than "paper" style, i.e. even rather trivial steps are written down 
explicitely. The goal is to pave a broad and smooth avenue to a deep under- 
standing of the smooth particle hydrodynamics method rather than just to 
provide a bumpy trail. Due to this pedagogical scope, the focus of this re- 
view needs to be clear-cut, but rather narrow: it only discusses the numerical 
solution of the inviscid hydrodynamics equations, in the Newtonian, special- 
rclativistic and general-relativistic (fixed metric) case. 

For practical astrophysical simulations often sophisticated additional physics 
modules are required and their implementation may pose additional numeri- 
cal challenges which arc beyond the scope of this review. In the following we 
provide a brief list of additional physics that has been implemented into SPH 
and point to references that are intended as starting points for further reading. 

• Gravity: 

Self-gravity is for many astrophysical problems a key ingredient. A straight 
forward calculation of pairwise gravitational forces between N particles re- 
quires a prohibitively large 0{N'^) number of operations and therefore such 
an approach is not feasible for realistic problems. The most natural gravi- 
tational force evaluation for a purely meshfree method like SPH is the use 
of a tree, either in the form of an oct-tree [7,3,8-12] or a binary tree [IS- 
IS] both of which require only 0(A^log(A^)) operations. Other possibilities 
include particle-mesh methods [16-18] in which the particles are mapped 
onto a mesh. The latter is used to efficiently solve Poisson's equation and to 
subsequently calculate the accelerations at the particle positions. Dehnen 
[19,20] has used ideas from the Fast Multipole Method (FMM) [21] for a 
fast tree code which scales proportional to 0{N). Several hybrid approaches 
which combine elements of different methods exist [22-26]. 
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• Equations of state: 

For some problems polytropic equations of state (EOS), p = Kp^, with K 
being a function of entropy, p the matter density and T the adiabatic expo- 
nent, are sufficiently accurate. Other problems require more sophisticated 
approaches. This is particularly true in cases where the matter compressibil- 
ity depends on the density regime, as, for example, in neutron stars where 
the effective polytropic exponent above nuclear matter density is F 2.5, 
but drops to values close to 4/3 at lower densities. Different nuclear EOSs 
have been used [27-31] in the context of neutron stars, at lower densities 
and in cases where the matter composition is important for other processes 
further physical EOSs have been applied [32-36] . 

• Solid state mechanics 

In the context of planetary impacts additional concepts such as material 
stresses, fracture physics or plasticity criteria need to be implemented into 
SPH. This branch of development was initiated by the work of [37] and fur- 
ther continued by [38-45] . More recently, also porosity models were included 
in SPH [46,47]. 

• Physical viscosity: 

For some applications the solution of the non-dissipative hydrodynamics 
equations is inadequate and physical (as opposed to artificial) viscosity 
needs to be modelled. Examples include the angular momentum transport 
in accretion disks, see e.g. [48], or the dissipation of sound waves as a mech- 
anism to heat intra-cluster media [49]. Various authors have implemented 
the viscous stress tensor into SPH [50-56]. This requires particular care in 
the implementation of higher order derivatives which -in a straight forward 
kernel approximation- can be sensitive to particle disorder. 

• Thermal conduction: 

Thermal conduction has, for example, been proposed as a possible heating 
mechanism for offsetting central cooling losses in rich clusters of galaxies. 
Its implementation poses some challenges due to second derivatives that 
enter the conductive term of the energy equation. Usually, a particluar dis- 
cretization of the second derivative due to [57] is used in modelling thermal 
conduction [58,59]. 

• Nuclear burning: 

Often the energy release due to nuclear reactions is too slow to influence gas 
flows substantially on a dynamical time scale. Under special circumstances, 
however, explosive nuclear burning can be triggered and in such cases the 
nuclear energy release can become the main driver of the dynamical gas 
evolution. In such cases the coupling between the hydrodynamics and the 
nuclear reactions is challenging due the huge difference in their intrinsic 
time scales. The flrst implementation of a nuclear reaction network into 
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SPH goes back to [60]. The 14 isotope a-chain network of this first study 
has also been used with updated nuclear reaction rates in studies of white 
dwarf coalescences [32]. More recently, further small networks [61,62] have 
been coupled with SPH [36] to study white dwarf coalescences [35,63], col- 
lisions between them [64,65] and tidal disruptions of white dwarfs by black 
holes [66,67]. 

• Chemistry: 

Chemical abundances carry the imprints of the past evolutionary processes 
that shaped today's galaxies. Their chemical enrichment is inseparably in- 
terweaved with the formation of structure in galaxies since -being closely 
related to supernovae- it drives local energy feedback and gas expansion, but 
on the other hand it also triggers -via metallicity-dependent cooling- the 
contraction and collapse of local gas structures. A large variety of chemical 
evolution approaches has been implemented into SPH, for a starting point 
one may want to consult [68-79] . 

• Radiation: photons and neutrinos 

Closely related to chemical species is the interaction between photons and 
the ambient gas. For some problems relatively simple cooling and heating 
recipes can be applied, e.g. [3,80,12]. Optically thick diffusion has in fact 
already been implemented in one of the very first SPH-papers [1]. Later, 
more sophisticated approaches to calculate the required second derivatives 
were suggested [57]. For a recent implicit diffusion approach, see [81]. Flux- 
limited diffusion schemes have been suggested by [82-84] , Monte-Carlo-type 
approaches have been followed by [85-88] and ray-tracing has been applied 
by [89-91]. Several further methods to treat radiation in SPH do exist [92- 
97]. 

When matter densities are so large that photons are "trapped" , neutrinos 
can become the only cooling agents. For white dwarfs and less compact 
stars matter is transparent to neutrinos and they can just be treated as a 
local drain of thermal energy. For higher densities, say in a neutron star or 
the inner regions of the disk around a hyper-accreting black hole, opacity 
effects may prohibit the free escape of neutrinos. A further example of the 
astrophysical importance of neutrino physics is the delayed explosion mech- 
anism [98,99] of core-collapse supernovae where the energy deposition via 
neutrinos is thought to be vital to re-energize the stalled shock after it has 
crossed the outer layers of the stellar iron core. In this context, multi-fiavour 
neutrino cooling and heating has been implemented into SPH [27,31] in an 
entirely particle-based, flux-limited diffusion approach. To model compact 
object mergers, neutrinos have been implemented into SPH via a hybrid, 
particlc-mcsh approach [100] and by using local disk scale heights to esti- 
mate neutrino opacities [33,34]. 
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• Magnetic fields: 

Magnetic fields are prevalent everywere in the cosmos, but they pose a noto- 
rious numerical challenge not least because of the difficulty in fulfilling the 

— * 

W • B — 0-constraint. The most obvious approach is a straight-forward SPH 
discretization of the ideal MHD equations [2,101-106]. While it performs 

— * 

reasonably well in some problems, particle disorder and finite V • B values 
can trigger numerical instabilities [105]. Several approaches have been devel- 
oped that apply regularization techiques to overcome the problems related 
to particle disorder [107-111]. Recently, an approach has been implemented 
[112,113] that advects so-called Euler-potentials, a and /3, with the SPH 
particles. The magnetic field can be reconstructed from the values of these 
potentials at the particle positions via B = V« x V/?. This approach guar- 
antees the fulfillment of the V ■ B = 0-constraint by construction and has 
shown excellent results in a large number of test cases [112]. 

• Further "subgrid" models: 

Often, numerical models of astrophysical processes require the feedback from 
physical processes on sub-resolution scales. Such processes are usually imple- 
mented as (parametrized) "subgrid" models that transmit the main effects 
from the unresolved processes to the fluid. Cosmic rays present one such 
example. They contribute significantly to the pressure of the interstellar 
medium in our Galaxy and therefore they may be important in regulating 
star formation during the formation and evolution of galaxies. For the im- 
plementation of the effects of cosmic rays into SPH see [114-118]. Another 
example is the modeUing of the feedback of (unresolvable) accretion-driven 
outfiows during the merger of galaxies to study the self-regulation of black 
hole growth [119,120]. In simulations on galactic scales stars -their birth, 
feedback via outfiows and possible supernova explosions- also need to be 
modelled as effective sub-scale models. By their very nature such approaches 
are closely linked to those of radiative transfer and chemical enrichment. For 
the various approaches we refer to the literature [121-143,12,144,145]. 

The remainder of this review is organized as follows. In Sec. 2 the most basic 
form of SPH ("vanilla ice") is derived by discretizing the equations of La- 
grangian hydrodynamics. This form is somewhat ad hoc, yet it performs well 
in practice and it is at the heart of many SPH codes that are regularly used 
in the astrophysics community. In Sec. 3, the Lagrangian of an ideal fiuid is 
discretized and subsequently used to derive a modern version of SPH. Since 
it follows from a variational principle, this formulation avoids any ambiguity 
with respect to a symmetrization of the equations and it naturally leads to 
additional terms that result from derivatives of the local hydrodynamic resolu- 
tion length, the so-called "grad-h-terms" , [146,147]. Armed with this gadgetery 
suitable Lagrangians are used to elegantly derive the SPH-equations both for 
the special- and the general-relativistic case, see Sec. 4. The summary and a 
brief outlook are presented in Sec. 5. 
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2 Basic concepts of Smooth Pcirticle Hydrodynamics 



2.1 Lagrangian hydrodynamics 



In contrast to grid-based (Eulerian) methods, Smooth Particle Hydrodynamics 
(SPH) is purely Lagrangian. In the Eulerian picture, derivatives are calculated 
at a fixed point in space while, in the Lagrangian description, they are evalu- 
ated in a coordinate system attached to a moving fluid element. Thus in the 
Lagrangian approach we follow individual fish rather that staring at the pond 
and watch the swarm pass by. The Lagrangian (or substantial) time derivative, 
d/dt, is related to the Eulerian time derivative, d/dt, by 

d dx^ d d _ ^ d 

— = \ = V • V H . (1) 

dt dt dx^ dt dt ^ ' 



Applied to the Eulerian continuity equation, 

f^+V-{pv) = 0, (2) 

one finds, using dp/dt = v ■ Vp -|- dp/dt, the Lagrangian form 

|.-PV... (3) 



The momentum conservation equation for a non-viscous fiuid, the so-called 
Euler equation, reads in Lagrangian form 



i.e. apart from "body forces" such as gravity or magnetic fields embodied in the 
quantity /, the fiuid is accelerated by gradients of the pressure P. The energy 
equation follows directly from the (adiabatic) first law of thermodynamics, 
Eq. (7), together with Eq. (3): 

du P dp Pyj ^ /rX 

— = — — = v-v. (5) 

dt p^ dt p ^ ^ 
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First law of thermodynamics 



For SPH one needs the first law of tliermodynamics, dU = dQ — PdV, in terms 
of "specific" quantities. Restriction to adiabatic processes allows to drop the 
dQ-term. On a "per mass basis" , the energy U becomes u ( "energy per mass" ) 
and the volume, V, becomes "volume per mass", i.e. with dV d{^/ p) — 
—dpjfP'. For the case without entropy generation the first law reads 

P 

du = -i^dp, (6) 
P 

from which 



follow. In a relativistic context it can be advantageous to work with quantities 
"per baryon" : 



where n is the baryon number density in the local rest frame. 



The set of equations (3), (4), (5) must be closed by an equation of state that 
relates quantities such as the pressure, P, or the speed of sound, c^, to macro- 
scopic fluid quantities such as density or temperature. All the microphysics, 
say whether the pressure is produced by collisions in a Maxwell-Boltzmann 
gas or by degenerate electrons due to the Pauli principle, is embodied in the 
equation of state. It can be as simple as a polytrope, or, for more complicated 
cases, say for hot nuclear matter, it may only be available in tabular form, see 
the introduction for a link to the existing literature. 

2.2 The SPH kernel interpolation 

In the following, discrete representations of the continuous Lagrangian hydro- 
dynamics equations are derived. In SPH, the interpolation points ( "particles" ) 
are moved with the local fluid velocity ^ , derivatives are calculated via a kernel 
approximation without the need for finite differences. In this way the partial 
differential equations of Lagrangian fluid dynamics are transformed into ordi- 



Variants that use local averages of fluid velocities exist, see the "XSPH" -approach 



du P dp 
dt p^ dt 



and 




(7) 




(8) 



of [148]. 
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nary differential equations. 

In principle, there is some freedom in discrctizing the fluid equations. But 
to ensure that the physically conserved quantities are also conserved by con- 
struction in the discretized particle equations the latter need to possess the 
correct symmetries in the particle indices. In the simplest, "vanilla ice" ver- 
sion of SPH, the symmetrization is imposed somewhat ad hoc, yet, it works 
well in practice and it is commonly used throughout astrophysics, see e.g. 
[4,149,150,135,14,30,151]. For different symmetrization possibilities see [3,152,50,153] 
and Sec. 2.5. 

2.2.1 Interpolating function values 

At the heart of SPH is a kernel approximation in which a function /(r) is 
approximated by 

= / f{f')W{f- h) dV, (9) 



where W is the so-called smoothing kernel (or window function) and the 

smoothing length, h, determines the width of this kernel. Obviously, one would 
like to recover the original function in the limit of an infinitely small smoothing 
region and therefore the kernel should fulfill 

lim fh{r) = /(r) and / W{r- P, h) = 1, (10) 
/i— »o J 



i.e. apart from being normalized, the kernel should have the ^-distribution 

property in the limit of vanishing smoothing length. 

To arrive at a discrete approximation, one can write the integral as 

/,(r) - / l^W{T-?,h) p{?)dh\ 
J p{r) 



(11) 



where p is the mass density and subsequently one can replace the integral by a 
sum over a set of interpolation points ("particles"), whose masses, nib, result 
from the p{r') d^r' term^ 

m^Y.—fbW(r-^>^^h). (12) 

b Pb 

^ For simplicity, wc arc omitting the subscript h in what follows. We also drop 
the distinction between the function to be interpolated and the interpolant, i.e. we 
use the same symbol / on both sides of the equal signs. Throughout the article 
subscripts such as Af, refer to the values at a particle position, i.e. Ai, = A{f},). 
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The summation interpolant Eq. (12) of the density then reads 

p{r^^^mbW{r-ri„h). (13) 

b 

This density estimate by summing up kernel-weighted masses in the neighbor- 
hood of the point f plays a central role in the derivation of the SPH equations 
from a Lagrangian, see Sec. 3. Note that h has not been specified yet, this will 
be addressed at a later point. In practice, the density can also be calculated 
via integration of a discretized version of Eq. (3) via Eq. (31), see below. In 
most applications the difference between integrated and summed density es- 
timates is completely negligible. The summation form, however, is certainly 
numerically more robust, integration can, in extremely underresolved regions, 
produce negative density values. The summation form does not assume the 
density to be a differentiable function while Eqs. (3)/(31) do so [154]. The 
standard practice in SPH is to keep the particle masses fixed so that the mass 
conservation is perfect and there is no need to solve the continuity equation. 



2.2.2 Approximating derivatives 

To calculate derivatives, one takes the analytical expression of the summation 
approximation, Eq. (12), 

V/(r-) = E —fb^W{r- n, h), (14) 

b Pb 



i.e. the exact derivative of the approximated function is used. The kernel 
function is known analytically, therefore there is no need for finite difference 
approximations . 

Processes such as diffusion or thermal conduction require second derivatives. 
One could now proceed in a straight forward manner and take another deriva- 
tive of Eq. (14). Yet, this estimate is rather sensitive to particle disorder and 
therefore not recommended for practical use. A better prescription is [57] 

(VV)a = 2E-(/a-A) — , (15) 



where Wab is related to the kernel by VWab = Gab Wab, Gab is the unit vector from 
particle b to particle a, iab = Tab/ Tab and fab = fa — H- For more information 
on higher order derivatives we refer to the existing literature [57,58,155,6]. 
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2.2.3 The kernel function 



To restrict the number of contributing particles in the sum of Eq. (12) to a lo- 
cal subset, the kernel should have compact support, otherwise the summation 
would extend over all n particles and produce a numerically intractable, com- 
pletely inefficient n^-method. Although for geometries such as flattened disks 
non-radial kernels seem like a natural choice [156-158] they have the disad- 
vantage that it is difficult to ensure exact angular momentum conservation, 
see below. Therefore, usually radial kernels with W{r — f^, h) = W{\r — f^\, h), 
are used in SPH. 

The accuracy of the kernel interpolation is in practice rather difficult to quan- 
tify unless strongly simplifying assumptions about the particle distribution 
arc made, which are usually not met in reality. One can expand f{r') in the 
integral approximant, Eq. (9), into a Taylor series around r 



By explicitly writing down the integrals over the different terms in the Taylor 
expansion on the RHS and comparing to the LHS, one can determine how 
well the approximation agrees with the original function. By requiring that 
error terms vanish, one can construct kernels of the desired order. In practice, 
however, such kernels may have unwanted properties such as negative values in 
certain regions and this may, in extreme cases, lead to unphysical (negative!) 
density estimates from Eq. (13). Although substantial work has been invested 
into constructing more accurate kernels [159,160,105,161], none has been found 
to overall perform substantially better in practice than the "standard" cubic 
spline SPH kernel [5] which is used in almost all SPH simulations. Particular 
higher-order kernels, however, may have advantages in the context of resolving 
Kelvin-Helmholtz instabilities near large density gradients [162]. In 3D, the 
cubic spline kernel reads 

1 - |g2 + |g3 for < g < 1 

1(2 -qf forl<g<2 , (17) 

for g > 2 



where q = \ f— r'\/h. This kernel is radial., i.e. it depends only on the absolute 
value of r — f*. With this kernel the integral approximation, Eq. (9), is 

A(r)=/(r-') + C/i2 + 0(/i^), (18) 
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where the qTiantity C contains the second derivatives of the function /. There- 
fore, constant and hnear functions are reproduced exactly by the integral rep- 
resentation Eq. (9). 

In practice, two approximations were applied: i) the integral interpolation, 
Eq. (9) and ii) the discretization, Eq. (12). The accuracy of the discretized 
equations depends on the distribution of the interpolation points ( "particles" ). 
In early work, error estimates were based on a purely statistical description 
assuming a random distribution of particles. In a simulation, however, the 
particle distribution is not random, but depends on both the kernel used and 
the dynamics of the considered system. For this reason, numerical experiments 
turned out to be far more accurate than predicted by these early estimates 
[6]. 

In the following box we collect some expressions for the derivatives of radial 
kernels that are frequently used throughout this text. 



2.3 The "vanilla ice" SPH equations 



2.3.1 Momentum equation 

One could now proceed according to "discretize and hope for the best" and 
this was historically the first approach. A straightforward discretization of 
Eq. (4) (no body forces) yields 

'| = -lj:!55nV„H/.. (19) 

dt Pa Y Pb 



This form solves the Euler equation to the order of the method, but it does 
not conserve momentum. To see this, consider the force that particle h exerts 
on particle a, 

- ( dva\ mamb 

Fba = ^a-jT = Pb^aWab, 20 

V dt J Pa Pb 



and similarly from a on 6 

Fab = rUb— = Pa'^bWba = Pa^ aWab, 21 

V dt J ^ Pb Pa Pa Pb 



where Eq. (26) was used. Since in general Pa 7^ Pb, this momentum equation 
does not fulfill Newton's third law ("actio = reactio") by construction and 
therefore total momentum is not conserved. 
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Derivatives of radial kernels 

Throughout the text, we use the notation fi,k = H — fk, r^fc = \ fhk\ and Vi,k = Vb — Vk- Derivatives 
resulting from the smoothing lengths are ignored for the moment, they will receive particular 
attention in Sec. 3. By straight- forward component- wise differentiation one finds 

or a \rb-rk\ 



(23) 



where ibk is the unit vector from particle k to particle b, 

d 1 _ ebkjSba - 5ka) 

o — * I — * — * I I — * — * I o 

ora\rb-rk\ In - Tk^ 

and 6ij is the usual Kronecker symbol. Another frequently needed expression is 

drab _ drab dxg ^ drab dya ^ drgb dza ^ drgb dxh ^ drab dijb ^ drgb dzb 
dt dxg dt dyg dt dzg dt dxb dt dyh dt dzb dt 

= ^aTab ■ Va + ^bTab ' Vb = ^ aTab ' Va " ^aTab ' Vb = ^ aTab ' Vab = ^gb ' Vgb, (24) 

where dvgb/dxb — —drab/dxg etc. was used. For kernels that only depend on the magnitude of 
the separation, W{fb — fk) — W{\fb — Htl) = Wbk, the derivative with respect to the coordinate 
of an arbitrary particle a is 

V7 u/ ^ u/ dWbkdrbk dWbk ^ 

^aWbk = -^Wbk = ^ = ^ ebkiha - ha) = ^bWkb{dba " ha): (25) 

org (JTbk OTg Orbk 

where Eq. (22) was used. This yields in particular the important property 

V7 ^ Ti/ ^^^b drgb dWgb . dWgb drab d 

VaW^ab ^ W^Wab= ^ = Cgb = — = - Wgb ^ -^bWgb- 26) 

drg drab dra drgb drgb drb drb 



The time derivative of the kernel is given by 

dWgb dWgb drgb dWgb (fa - fb) ■ ( fa - Vb) dWab 



dt drab dt drgb rgb dr. 



ab 



eab ■ fab = fab ' a^gb- (27) 



But a slightly more sophisticated approach yields built-in conservation. If one 
starts from 



solves for VP/ p and applies the gradient formula (14), the momentum equa- 
tion reads 
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dtpab b Pb Pb 

— Emb(^ + ^)V.W^. (29) 

b \Pa Pb/ 

Because the pressure part of the equation is now manifestly symmetric in 
a and b and VaWab = — V^Wfea, the forces are now equal and opposite and 
therefore total momentum and angular momentum (see below) are conserved 
by construction. It is worth pointing out some issues: 

• So far, any complications resulting from variable smoothing lengths have 
been ignored. Conservation is guaranteed if a) constant smoothing lengths 
are used (which is a bad idea in practice!) or b) "VaVFai" is symmetric in the 
smoothing lengths, which can be achieved, for example, by using hab = {ha + 
hb)/2 [4] or by replacing the gradient via \SoW{rab,K) + ^ aW{rab,hb)\/2 
[3]. 

• The conservation relies on the force being a term symmetric in a and h times 
Cab (keep in mind that V aWah oc Cab, sec Eq. (26)). 

• SPH's success largely relies on its excellent conservation properties which is 
guaranteed by the correct symmetry in the particle indices. 



2.3.2 Energy equation 

A suitable energy equation can be constructed from Eq. (5) in a straight 
forward manner: 

dUa Pa dpa Pa d ( ^ \ Pa ^ ^ ^ 

2^ rribWab = ^ 2^ mbVab ■ Va^Kfe, (30) 



dt pI dt pI dt Vt / pI 



b 



where Eq. (27) was used. Together with an equation of state, the equations 

(13), (29) and (30) form a complete set of SPH equations. 

For later use we also note that the velocity divergence can be conveniently 

expressed by using Eqs. (3) and (27) as 

(V • i;)„ = - -^ = -1 ^m,i;„, • V„W^„,. (31) 

Pa dt Pa Y 



One can derive an alternative energy equation, and this will allow a smoother 
transition to the rclativistic equations, by using the specific "thermokinetic" 
energy ia = Ua + \v1 instead of the specific thermal energy The corre- 
sponding continuous evolution equation, 

§^-lv.(PS), (32) 
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can be written as 



| = _^V.(H-.-.v(g. (33) 



Applying Eq. (14) yields 



dt PabPb b Pb Pb 

= - E f ^ + ^ V ^aWab. (34) 
b \ Pa Pb J 

As we will see later, see Eq. (167), this is similar to the relativistic form of the 
energy equation. 



2.4 Conservation properties 



We can easily check the numerical conservation of the physically conserved 
quantities. The change of the total particle momentum is 



i:rr^a^=j:j:^ba=Uj:j:^ba+j:j:^ba] 

a a b ^ \ a b a b / 

= \\Y.Pab + Y. Fba] = \ ( E(^«f + Fba) \ = (35) 
\ a,b b,a ) \ a,b 

where the (dummy) summation indices were relabeled after the third equal 

sign and the notation from Eq. (20) and Fab = —Fba were used. 

The proof for angular momentum is analogous. The torque on particle a is 

Ma^faXFa^faX (tuJ-^^ = X (36) 



and thus 



^ = J2Ma = E^« X Fba=^ (E^« X + X Fba 

"'t a a,b ^ \a,b a,b 

= I lE^"" X ^ba + E^b XFab]^^ [Y^ifa - fb) X | = 0. (37) 

\ a,b b,a I \ a,b 
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Again, the dummy indices were relabeled and Fab — —F^a was used. The 
expression finally vanishes, because the forces between particles act along the 
line joining them: Fab oc VaWab oc iab oc {fa - n), see Eq. (26). 
The total energy changes according to 



d -r^ f ,1 V- d^a ^ fdua dva\ , . 

dt^\ 2 ''J ^ " dt ^ \dt dt I ^ ' 

Inserting Eqs. (30) and (29) we have 



dE 
~dt 



dE ^ 
— — = y nia 
dt V 



^JZ'^bVab • ^aWab " " "^M "| + "l ) ^ a^ab 
Pa b b \Pa PbJ 



Pa Pa 
= Y.'^amb-^Va ■ "^aWab - Xl"^a"^6^V6 • V aWab 

a,b Pa a,b Pa 

P Pb 
-^malTlb-^Va ■ VaWab " Y.'^o-'^b^^a ' Val^ab 

a,b Pa a,b Pb 

= - E marub f ^ + ^ V WaWab. (39) 

\ pI Pb J 

Prom this equation we could have directly read off the evolution equation for 
the specific thermokinetic energy, Eq. (34). To show the conservation of the 
total energy we can apply the same procedure as in the case of the momentum 
equation since we have again a double sum over a quantity symmetric in the 
indices a, b multiplied by a quantity which is anti-symmetric in the indices 
a, b. Therefore, the double sum vanishes and total energy is conserved. 



2.5 Alternative SPH discretizations 



Like other numerical schemes, SPH offers some freedom in discretizing the con- 
tinuum equations. At this point it is worth remembering that we had enforced 
momentum conservation by using a particular form of the gradient on the 
RHS of the momentum equation, see Eq. (28). Following the same reasoning 
one could as well have used [5] 



VP 
P 

which yields 

dv. 




(40) 




VaWab (41) 
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for the momentum equation. This form is symmetric in the particle indices for 
any value of A and therefore also conserves momentum by construction. This 
equation is consistent with the energy equation 



dUa Pa ^ Vab ^ „ , 
Pa b Pb 



(42) 



and the continuity equation 

dpa 2-A V7 tj/ 

Pa Z^mb^^-VaWab, 



dt 



Pb 



(43) 



which is the SPH-discrete form of 
dpa 2-X 



dt 



P 



(44) 



The consistency of the above SPH equation set can be shown by using a 
generalized variational principle [105] that does not assume that the density 
is given as a function of the coordinates, but instead assumes a given form of 
the continuity equation. 

If A = 1, the continuity equation only depends on the particle volumes rUb/ Pb 
rather than the particle masses. This can reduce numerical noise in regions 
of large density contrasts where particles of very different masses interact 
[163]. Marri and White [137] found improved results in their galaxy formation 
calculations by using (among other modifications) a value of A = 3/2. 
The SPH equations can be generahzed even further after noting that the 
continuity equation can be written as [105] 

'*'' = *L-.Vm-V.fa|, (45) 



dt V* 



with ^! being a scalar variable defined on the particle field. This leads to SPH 
equations of the form 
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The "vanilla ice" equation set, Eqs. (29) and (30), is recovered for ^ = 1, 
choosing ^ = pj^fP one recovers the momentum equation of Hernquist and 
Katz [3]. 

The above equation set has been used in recent work [162] to address SPH's 
weakness of resolving Kelvin- Helmholtz instabilities across large density jumps. 
The authors choose a different auxiliary function for each of the above equa- 
tions, and perform an error and stabihty analysis leading them to 
suggest = ^'^ = and = p. Here A is the entropy-dependent part of 
the polytropic equation of state, P — A{s)p^ and F the polytropic exponent. 
This prescription leads to a sharper density transition that is more consistent 
with the entropy transitions between both parts of the fluid. Together with a 
higher-order kernel and large neighbor numbers they find convincing results 
in Kelvin-Helmholtz instability simulations. 



2.6 Adaptive resolution 



Whenever densities and length scales vary by large amounts, the smoothing 
length h should be adapted in space and time. Several ways to adjust the 
smoothing lengths have been used over the years [2,164,165,3,4,166]. One of 
them seeks to keep the number of neighbors of each particle (approximately) 
constant [3] , another to integrate an additional equation that makes use of the 
continuity equation [4] . One can use the ansatz 

m^(l!LX" (49) 

where the index labels the quantities at the beginning of the simulation. 
Taking Lagrangian time derivatives of both sides together with Eq. (3) yields 

^ = lhaiV . V)a. (50) 

Another convenient way is to evolve h according to 

ha = v( — Y\ (51) 



where r) should be chosen in the range between 1.2 and 1.5 [165,105,6]. 

Note, that the above SPH equations were derived under the assumption that 
the smoothing lengths are constant. Varying the smoothing lengths while using 
the above equation set is strictly speaking inconsistent. A consistent formu- 
lation that involves extra-terms will be discussed in detail in Sec. 3.3. The 
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importance of these corrective, or "grad-h" terms is problem- and resolution- 
dependent [146,112]. 

2.1 Artificial dissipation 

2.7.1 Purpose and general reasoning 

In gas dynamics, even perfectly smooth initial conditions can steepen into 
discontinuous solutions, or "shocks", see e.g. [167-169], which are nearly om- 
nipresent in astrophysics. On length scales comparable to the gas mean free 
path these solutions are smooth (as a result of the physical viscosity that is 
always present to some extent), but on the macroscopic scales of a simulation 
(which is usually orders of magnitude larger) , the very steep gradients appear 
as discontinuous. 

The strategies to deal numerically with shocks are (broadly speaking) two-fold: 
i) one can either make use of the analytical solution of a Riemann problem 
between two adjacent computational entities (either cells or particles) or ii) 
broaden the discontinuity to a numerically resolvable length across which gra- 
dients can be calculated. The latter is instantiated by adding extra, "artificial" 
dissipation to the flow. This can be achieved by choosing particular numerical 
discretization schemes that in fact solve continuum equations that are only 
similar to the original ones, but contain additional, higher-order deviatoric 
terms, such as, for example, the Lax scheme [170]. Alternatively, one may ex- 
plicitely add prcssurc-like terms to the fluid equations. The addition of such 
"ad-hoc" terms is one of the oldest techniques [171] in the relatively young field 
of computational physics. Most astrophysical SPH implementations follow the 
latter strategy, e.g. [3,14,12,112], but variants that use Riemann-solvers do also 
exist [172,173]. 

Strictly speaking, real discontinuities with their infinite derivatives arc no 
proper solutions of the ideal hydrodynamics equations Eqs. (3)- (5). Artifi- 
cial viscosity (AV) aims at spreading such discontinuities over a numerically 
resolvable length. In the words of von Neumann and Richtmyer (from their 
seminal paper [171]), the "idea is to introduce (artificial) dissipative terms 
into the equations so as to give the shocks a thickness comparable to (but 
preferentially larger than) the spacing ... [of the grid points]. Then the dif- 
ferential equations (more accurately, the corresponding difference equations) 
may be used for the entire calculation, just as though there were no shocks 
at all." Although guidance by physical viscosity can sometimes be helpful, 
artificial viscosity is not meant to mimic physical viscosity, it is instead an ad 
hoc method to produce on a resolvable scale what is the result of unrcsolvable, 
small-scale effects. In this sense one can think of it as a kind of sub-grid model. 
AV should have a number of desirable properties [174] and must not introduce 
unphysical artifacts. It should always be dissipative, i.e. transfer kinetic en- 
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crgy into internal one and never vice versa, which is not a completely trivial 
task in 3D. AV should be absent in rigid and (shockless) differential rotation 
and in uniform compression. In other words, it should be "intelligent" in the 
sense that it distinguishes uniform compression from a shock. This can be 
done either via a tensor formulation of artificial viscosity, e.g. [175], or via 
"limiters" that detect such motion and subsequently suppress the action of 
AV. Generally, AV should go smoothly to zero as the compression vanishes 
and it should be absent for expansion. Moreover, to not destroy one of SPH's 
most salient strength, it must be implemented consistently into momentum 
and energy equation, so that the conservation of energy, momentum and an- 
gular momentum is still guaranteed. 

2.1.2 Bulk and von-Neumann-Richtmyer viscosity 

Guided by the requirements that i) no real discontinuities occur, ii) the thick- 
ness of the "shock layer" is everywhere of the order of the resolvable length 
scale iii) no noticeable effects occur away from shock layers and iv) the 
Rankine-Hugoniot relations [167] hold when considering length scales that are 
large in comparison to the thickness of the shock layers, von Neumann and 
Richtmyer suggested an artificial pressure of the form 

qNR^C2pf{V-v)^ (52) 

where C2 is a dimensionless parameter of order unity, to be added to the 
hydrodynamic pressure. Their approach yielded reasonably good results in 
strong shocks, yet it still allowed oscillations in the post-shock region to occur. 
To avoid them, usually an additional term [176] is introduced that vanishes 
less rapidly and has the form of a bulk viscosity [167] 

?b = -Cipcs/(V • v). (53) 

Here Cg is the sound velocity and Ci is a parameter of order unity. Often, the 
two contributions are combined into an artificial viscous pressure, 

gvisc = -cipcJiW ■ v) + C2pf{V ■ vf. (54) 

The term I ( V-v) is thereby an estimate for the velocity jump between adjacent 
cells or particles. 

2.7.3 "Standard" SPH viscosity: reasoning and limitations 

A straight-forward approach to transform these ideas into an SPH artificial 
viscosity prescription would be to insert Eq. (31) into (54), to use particle 
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properties for the densities and sound velocities and identify the resolution 
length scale / with the smoothing length. Similar approaches have indeed been 
pursued [3], but we will focus here on an approach due to [177] that is proba- 
bly the most wide-spread form of SPH artificial viscosity and therefore often 
referred to as the "standard viscosity" . 

Instead of just increasing the hydrodynamic pressures, Pa, by viscous contri- 
butions, gvisc, one augments the pressure terms in the momentum equation 
(29) by an artificial contribution Uab- 




For simplicity, let us consider the bulk viscosity contribution first and re- 
strict ourselves to ID. The bulk contribution to Uab is then of the form 
—ciCsih/ p) dv/dx and a Taylor expansion of the velocity field, v{xa+{xb—Xa))i 
yields 




Xb - X. 



^ + O {{Xb - Xaf) . (56) 



Inserting this approximation, replacing particle properties by their averages, 
Aab = {Aa + Ab)/2, using the notation Xab = Xa — Xb and Vab — Va — Vb and 
replacing the possibly diverging denominator, 

— ^ 2 -u2 ^ (57) 
Xab Xg^fj -\- erig^f, 



the SPH prescription for the bulk viscosity reads 

-Ci^/Iab for Xab Vab < 



ria&.bulk — < 



Pab 





otherwise 



, where 



hab -^ah ^ab 



ab 



xlb + e/i,^ 



(58) 



ab 



The product Xab Vab thereby detects whether particles are approaching (< 0) 
and only in this case artificial viscosity becomes active. The quantity Hab 
has taken over the role of the term 1{S/ ■ v) in the original prescription of 
von Neumann and Richtmyer. Having said this, we can now treat the von 
Neumann-Richtmyer term, g^R, in exactly the same way. Adapting to the 
usual SPH notation, ci — > a, C2 — > /3, the artificial viscosity term reads 



H, 



ab 



H, 



afe.bulk 



H, 



a6,NR 



Pab 





for Tab - Vab <0 

otherwise 



(59) 
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with 




Fig. 1. Scliematic illustration of the shock tube problem. 



hab'f^ lb ■ Vab 
''"ab + ^^ab 



(60) 



Note that the scalar quantities have now been replaced by vector quantities 
for the use in 3D. Numerical experiments suggest the following values for the 
involved parameters: a ~ 1, /3 ~ 2 and e ~ 0.01. 
Accounting for artificial viscosity, the momentum equation reads 



dt 




+ Uab VaW, 



ab- 



(61) 



To have a consistent formulation, the energy equation must be modified ac- 
cording to 

dua Pa 1 

--jT = ^bVab ■ VaW^ab + ^ TTlbllabVab " VaM^ab- (62) 

dt Plk ^b 



This combination still conserves energy, linear and angular momentum by con- 
struction, the proofs are analogous to the ones outlined above. 
To illustrate the performance of this artificial viscosity prescription, let us 
consider a standard test case for hydrodynamic schemes, the so-called "Sod 
shocktube" [178] which is a particular realization of a Riemann problem. A 
tube filled with gas, such as illustrated in Fig. 1, contains a high density 
(pi = l)/high pressure {pi = 1) region that is separated by a wall from a 
low density (p2 = 0.25)/low pressure (p2 = 0.1795) region. For the Sod test a 
polytropic equation of state with exponent F = 1.4 is used. Once the wall is 
removed, the discontinuity decays into two elementary, non-linear waves that 
move in opposite directions: a shock moves into the unperturbed original low 
density region, while a rarefaction wave travels into the original high-density 
region. Between the shock and the tail of the rarefaction wave (the rightmost 
end of the rarefaction) two new states develop which are separated by a con- 
tact discontinuity across which the pressure is constant. The exact solution is 
sketched in Fig. 2. 

If no artificial viscosity is applied in such a problem, strong post-shock os- 
cillations occur as illustrated in Fig. 3 for the velocities in the Sod test. If 
the standard artificial viscosity prescription, Eq. (59) is used, the numerical 
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Fig. 2. The exact solution of the non-relativistic shock tube problem. 




Fig. 3. Velocity snapshot of the Sod shock tube problem. Shown is the exact solution 
together with a numerical solution without artificial viscosity. 

solution agrees very well with the exact one, see Fig. 4, provided that the 
interesting regions of the flow are adequately resolved ^ . 
While this artificial viscosity form performs very well in ID, it suffers from 
deficiencies if used in multi-D: i) for fixed parameters a and (3 it may affect 

^ Consistent with the finite resolution width of numerical shocks, we set up the 
initial density distribution as a Fermi function with a transition width of Ax: p{x) = 
{pi — 92)1 (l + 6xp(^^^)) + p2, where xs is the initial position of the shock. For 
the presented test Ax was set to 1.5 times the average particle separation. 
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density velocity 




-0.4 -0.2 0.2 0.4 -0.4 -0.2 0.2 0.4 




Fig. 4. Sod's shock tube: the exact solutions are shown by the red Hne, the numerical 
solution obtained with SPH are shown as circles. For this simulation 1000 particles 
were used, density was calculated via summation, and the smoothing lengths were 
updated according to ha = lA{ma/ Pa)- For this test a second-order Runge-Kutta 
integrator was used. 

the flow even if it is not really needed and ii) under certain conditions it can 
introduce spurious shear forces. Consider for example an idealized shear flow 
as sketched in Fig. 5. Assume that the velocity decreases vertically as sketched 
on the left and no shocks occur. For such a situation no artificial viscosity is 
needed. Nevertheless, as sketched for two example particles ("1" and "2") the 
scalar product r^b ■ Vah is finite and thus, via Hah in Eq. (60), introduces a 
viscous force that is unwanted. Similar configurations are encountered for the 
(astrophysically important) cases of accretion disks. 



2.7.4 Reducing artificial viscosity where unnecessary 

Several recipes were suggested to cure the deficiencies that the standard arti- 
ficial viscosity introduces in multi-D. 

One such recipe is the so-called "Balsara- switch" [179]. Its strategy is to de- 
fine a "limiter" that distinguishes shock from shear motion and suppresses AV 
in the latter case. One defines a quantity 

f ^ l(V-^)a| /goN 

\{V-v)a\ + \{Vxv)a\+O.OOOlCs,a/ha ^ ' 
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Fig. 5. Illustration how the "standard" SPH artificial viscosity introduces spurious 
shear forces: in this pure shear fiow the difference position vector has a finite pro- 
jection on the difference velocity vector (red) and thus introduces unwanted forces. 

and, as before, applies a symmetrized average of the limiters /, 11'^^ = Ilaf, fab- 
In pure compressional motion (|(V ■ w)| 7^ and |(V x -u)! = 0) the limiter 
reduces to unity and the standard viscosity is recovered, and in pure shear flow 
(|(V ■ ^7)1 = and \{V x v)\ 7^ 0) the action of artificial viscosity is suppressed 
{fab 1). This limiter has been found very useful in many cases [149,180,181], 
but it reaches its limitations if shocks occur in a shearing environment such 
as in an accretion disk [175]. 

Morris and Monaghan [53] suggested the use of time- dependent artificial vis- 
cosity parameters. The main idea is to have a and f3 only at non-negligible 
levels where they are really needed. They suggested to fix to 2a, to assign 
each particle its individual parameter, a^, and to evolve it according to an 
additional differential equation, 

dcXa Ola Q^min , ri /i^a\ 

-jT = + Sa, (64) 

dt Ta 



SO that Ua decays exponentially with an e-folding time Ta to a minimum value 
ainin, unless it is triggered to rise by the source term Sa- The time scale Ta 
should be chosen so that the viscosity persists for a few smoothing lengths 
behind the shock. This can be obtained by 



hn 



(65) 



where ^ ^ 0.1 ^ . Morris and Monaghan suggested Sa = max 



the source term. If one wishes to restrict the growth of a to Omax, one can 



-(v-t;)a,o 



for 



use Sa = max — (V 



.0 



[181]. Typical values for the numerical 



Note the similarity to the SPH-version of the Courant time step criterion. 
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parameters are amin = 0.1 and a^ax = 1-5. In a Sod test this prescription leads 
to non-negligible values of aa only in the vicinity of the shock front, elsewhere 
artificial viscosity is practically absent. This prescription has largely removed 
unwanted effects from AV in three-dimensional SPH simulations [181,182]. 
Note, however, that in homologous flow, v cxf, aa can still rise although it is 
not needed. This is counterbalanced to some extent by the increase of Cs,a ^-nd 
the decrease of ha, so that -via r^- the decay term in Eq. (64) becomes more 
dominant, but there is certainly still room for improvement. 



2. 7. 5 New forms inspired by Riemann solvers 

Monaghan used the analogy to Riemann solvers to motivate a new form of 
artificial dissipation which involves signal velocities and jumps in variables 
across characteristics [183]. The main idea of these "discontinuity capturing 
terms" is that for any conserved scalar variable A with Y^a'^^adA^/ dt = a 
dissipative term of the form 

f^l =Y.^i,^^^{Aa-A,)eab-VaWab (66) 



should be added, where the parameter aA,b determines the exact amount of 
dissipation and Vsig is the maximum signal velocity between particle a and b. 
Applied to velocity and thermokinetic energy this yields 



^ aVsig{Va - Vb) ■ e-ab^ 



diss 6 P"-b 
"a ~ H 



diss 6 



Pab 



eab ■ VaWaft, 



(67) 
(68) 



where, following [154], the energy including velocity components along the line 
of sight between particles a and 6, e* = \avsig{va ■ Cahf + o-u^sigUai has been 
used and different signal velocities and dissipation parameters were explicitely 
allowed for. If one uses dua/dt = de-a/dt — Va ■ dva/dt, one finds the dissipative 
terms of the thermal energy equation 



'dUg 

dt 



diss 



rrib 

Pab 



1 



aVsig -[Vab ■ Cab) + an<ig(Ma 



Ub) 



eab ■ VaW„6.(69) 



The first term in this equation bears similarities with the "standard" artificial 
viscosity prescription, see Eq. (62), the second one expresses the exchange of 
thermal energy between particles and therefore represents an artificial ther- 
mal conductivity which smoothes discontinuities in the specific energy. Such 
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artificial conductivity had been suggested earlier to cure the so-called "wall 
heating problem" [184] . Tests have shown that artificial conductivity substan- 
tially improves SPH's performance in simulating Sedov blast waves [112]. 
For non-relativistic hydrodynamics the maximum signal velocity between two 
particles can be estimated as [183] 

^^sig = Cs,a + Cs,b - Vab ' ^ab, (70) 

where Cs,fe is the sound velocity of particle k. Price [154] had reahzed that 
SPH's difficulty to treat Kelvin-Helmholtz instabilities across contact discon- 
tinuities with large density jumps [185] is closely related to a "blip" that occurs 
in the pressure at the contact discontinuities ^ . He suggested to use artificial 
conductivity only to eliminate spurious pressure gradients across contact dis- 
continuities and to this end suggested 



Clearly, this quantity has the dimensions of a velocity and vanishes in pressure 
equilibrium. This approach has substantially improved SPH's ability to treat 
Kelvin-Helmholtz instabilities [154]. 

To avoid conductivity where it is unwanted, one can follow again a strat- 
egy with time-dependent parameters. For the artificial viscosity one can use 
Eq. (64), and proceed in a similar way for One can use the second deriva- 
tive of the thermal energy, 

S„='^, (72) 

to control the growth «„. The second derivative can be calculated as in 
Eq. (15), the parameter e avoids that S^^a diverges as — > 0. 
According to a recent analysis [162], SPH's difficulty to treat Kelvin-Helmholtz 
instabilities results from a mismatch in the sharpness of pressure and density 
across the density jump. This can be either cured by generating entropy at 
the boundary and thus smoothing the pressure as in [154], or by obtaining a 
sharper density estimate. By a combination of using the freedom in discretiza- 
tion, see Sec. 2.5, a particular, higher-order kernel and an entropy-weighted 
density estimate together with large neighbor numbers, [162] also find con- 
vincing results Kelvin-Helmholtz instability simulations. 



^ It is not visible in our Sod shock tube in Fig. 4 since we had started from smoothed 
initial conditions. 




(71) 
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2.8 Time integration in SPH 



To integrate the ordinary differential equations of SPH one has to find a rea- 
sonable tradeoff between accuracy and efficiency of an integrator and most 
often the available computer power is better invested in larger particle num- 
bers rather than in high-order integration schemes. Since the evaluation of 
derivatives is usually very expensive, and in particular so for self-gravitating 
fluids, one tries to minimize the number of force evaluations per time step, 
which gives preference to low-order integrators. If the storage of derivatives 
(from earlier time steps) is not a concern, one can resort to multi-step methods 
such as Adams-Bashforth-type integrators, e.g. [170,186], giving higher-order 
time integration accuracy at moderate costs of force evaluations. 
After a short collection of commonly used time step criteria, we want to discuss 
briefly two methods that we find particularly useful: the Stormer-Verlet/leap 
frog algorithm, which appeals by its conservation properties and the class of 
Fehlberg methods whose advantage is the appropriate choice of the time step 
size based on monitoring the quality of the numerical solution. 



2.8.1 Time stepping 

We will briefiy collect here commonly used time step criteria, depending on the 
considered physical system, further criteria that capture the specific physical 
time scales have to be added. A criterion that triggers on accelerations, Oa, is 
[5] 

Ai/,a OC ^Jha/\aa\, (73) 



and 



^tcv,a OC \^ T (74) 



is a combination of a Courant-type ^ and viscous time step control, where 
is the quantity defined in Eq. (60). If one wishes to restrict the change of the 
smoothing length over a time step, one can use additionally[151] 

Ai,,„oci^. (75) 



® The Courant or Courant-Friedrichs-Levi (or CFL for short) criterion ensures that 
the numericai propagation speed of information does not exceed the physicai one. If 
a spatiai scafe of Ax can be resofved, the numericai time step has to be Ai < Ax/cg 
to ensure numerical stability[170]. 



29 



In a simulation, the minimum of the different time step criteria (with a suit- 
ably chosen prefactor) determines the hydrodynamical time step. If further 
physical processes, say nuclear burning, occur on much shorter time scales, 
one may resort to an "operator splitting" approach and integrate different 
processes with separate integration schemes, e.g. [67]. 

For inexpensive test problems or cases where the physical time scales are (more 
or less) the same for each particle, it is easiest to evolve all the particles on 
the same (shortest particle) time step. In many astrophysical examples, say 
cosmological structure formation, the collapse of a molecular cloud or the tidal 
disruption of a star by a black hole, the required time steps in different parts of 
the fluid may span many orders of magnitude. In such cases it is beneficial to 
group particles into block time steps of At„ = 2"Atinin, where n is an integer 
and Atmin is the smallest required time step, and to evolve each group of parti- 
cles separately. In this way the number of (expensive) force evaluations can be 
reduced by orders of magnitude. Successful implementations of such individual 
time step schemes can be found, for example, in [187,188,3,123,189,12,190,151]. 



2.8.2 Stdrmer-Verlet and leap frog 

The Verlet or leapfrog algorithm^ is particularly appealing due to its exact 
time reversibility. The original Verlet algorithm [191] is very easy to derive: 
start from two Taylor expansions for f{t+At) and r(t— At), add them and solve 
for f{t + At) to find the position update prescription of the Verlet algorithm 

f{t + At) = 2f{t) - f{t - At) + a{t)At'^ + 0{At^). (76) 



It is interesting to note that the position is updated without using the veloci- 
ties. But this comes at a price: for a position update two positions at earlier 
time steps are needed. The velocity can be reconstructed from the positions 
via centered finite differences: 

, f(t + At)-f(t-At) 
m = ^ ^ + 0{At'). (77) 



This simple integrator is only second-order accurate, but it is time reversible 
and has excellent conservation properties. To kick off a simulation at t = 0, 
one needs f{—At). It can be obtained by solving for f{—At) after inserting 
t = into Eq. (77): 

f{-At) = f{At) - 2Atvo, (78) 

''' This algorithm has many names: in astronomy it is often called Stdrm,er method^ 
in molecular dynamics it is usually named Verlet method in the context of partial 
differential equations it is usually referred to as leap-frog method. 
Note that the scheme needs modification if the acceleration depends on the velocity. 
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r(t + At) =f{t) + Atv[t + —). (81) 



where vo — v{t — 0). Inserting this into Eq. (76) provides the position update 
r{At) ^fo + Atvo+ ^ao At^, (79) 

which looks hke the first terms of a Taylor expansion around t = 0. 

The "leapfrog form" is obtained by defining velocities at half-steps, again via 

centered differences 

^/ At\ f(t)-f(t-At) , ^/ At\ f(t + At)-f(t) 

The last equation can be solved for the leapfrog form of the position equation 

At\ 

Starting from Eq. (76) and inserting Eq. (80) we find the velocity update of 
the leapfrog algorithm 

+ t) ^ " t) + ^^^^ 



So positions and accelerations are always evaluated at "full" time steps t", 
while the velocities are evaluated in between, either at i""^/^ or so the 

velocities are always "leaping" over the positions. If velocities at t are needed, 
they can be calculated according to 

v(t + f)+v(t-f) 
v{t) = ^ ' ^ ^ ^ ^. (83) 



The "Stormer-Verlet form" is particularly useful, since all quantities are eval- 
uated at the same point of time. Start from Eq. (77) and solve for 

f{t-At)=f{t + At)-2Atv{t). (84) 



Insert this into the Verlet position update, Eq. (76), to find the position update 
of the velocity Stormer-Verlet algorithm 

f{t + At) = f{t) + v{t) At + l-a{t)At'^. (85) 

To find the velocity update for the algorithm start from Eq. (77) and insert 
the Verlet position update, Eq. (76), 
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v{t) 



f{t) - f{t - At) 
At 



+ 



a{t) 



2 



At. 



(86) 



Completely analogous one finds for the velocity ai t + At 



v{t + At) 



fit + At) - fit) 

Xt 



+ 



dit + At) 

2 



At. 



(87) 



Adding the last two equations and using Eq. (77) gives 



fit + At) = f{t) + At 



d{t) + d{t + At) 
2 



(88) 



Together with Eq. (85) this update forms the Stormer-Verlet algorithm. It 
is a simple example of a "geometric integrator" which preserves qualitative 
features of the exact flow of the ODE by construction[192,193]. In particular, 
it conserves angular momentum exactly ^ . 

The integration of an elliptical Kepler orbit represents a simple, yet significant 
numerical experiment. As long as no deviations from purely Newtonian point 
mass gravity are considered, bound orbits are closed ellipses, e.g. [194,195], 
any "non-closure" is a numerical artifact due to finite integration accuracy. 
To define the problem, we use units in which a test body ( "planet" ) is accel- 
erated by 



and, since angular momentum is conserved in a central force field, one only 
needs to consider two spatial dimensions. We choose Tq = (1,0) and Vq — 
(0, 0.2) as initial conditions. We perform a simulation up to i = 150 with a 
slightly too large time step. At = 10"^. For comparison we also use a fourth- 
order Runge-Kutta scheme with the same time step. The results are displayed 
in Fig. 6. Due to the inappropriate time step the numerical inaccuracies in 
the velocity Stormer-Verlet case lead to a perihelion shift, but otherwise the 
behavior is physical: the test particle orbits the center of gravity on an eccentric 
orbit. In the Runge-Kutta case a qualitatively different and unphysical effect 
occurs: energy and angular momentum are dissipated for purely numerical 
reasons. The evolution of angular momentum is shown in Fig. 7. The "steps" 
in the Runge-Kutta case occur at perihelion since here the time steps are much 
too large. Therefore, the test particle drifts into a more circular, lower energy 
orbit. This illustrates that exact conservation can be more important that high 
order! 



For an elegant proof based on Newton's Principia sec [193]. 




r 



(89) 
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Fig. 6. Integration of the test particle orbit up to t = 150 for an eccentric or- 
bit {vo = (0,0.2), At = 10~^). Left: fourth-order Runge-Kutta. Right: velocity 
Stormer- Verlet . 



0.201 




Fig. 7. Angular momentum evolution of the long-term comparison. While the veloc- 
ity Stormer- Verlet method conserves angular momentum to machine precision, the 
fourth-order Runge-Kutta method loses about 2 percent. As a result of this purely 
numerical dissipation the planet spirals inward. 

2.8.3 Time stepping with quality control: Fehlberg methods 

The size of the time step At is crucial for the accuracy (and stability!) of the 
numerical solution of an ODE. It has to be a small fraction of the shortest 
physical time scale of the problem at hand. For a complex multi-physics sys- 
tem this may require to determine the time scales of every involved physical 
process. In the worst case, these may not be known and, on top, the meaning 
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of "a small fraction" may subject to trial and error. 

Ideally, a method should decide itself which step size to take. This is the phi- 
losophy behind adaptive step size control or Fehlberg methods [196-198,170]. 
The main idea is to use two estimates of different order for the solution and 
use them to estimate the error growth rate. Based on the comparison of the 
estimated and the tolerable growth rate, the chosen time step may be accepted 
or rejected. Independent of whether the time step was accepted or not, the 
growth rate can provide an educated guess for the size of the next time step. 
The procedure is the following: 

• Choose a trial time step At. 

• Use two different methods to obtain results at the end of the time step, ai 
and 02. To be reasonably efficient, the two methods should use nearly the 
same evaluations of the RHS of the ODE. 

• Use Oi and 02 to calculate an approximate local truncation error E. 

• Compare the error growth rate e ^ E/At with some acceptable error growth 
rate, e'^'^'^. If the rate is small enough, e < e'^'^'^, accept the time step, t""*"^ = 
t"' + At and i/""*"^ = 02, where 02 is the more accurate method. Otherwise, re- 
take the step with a reduced step size determined from e and the previously 
tried step size. 

Obviously, the above strategy can be applied to a large variety of integra- 
tor combinations. An often used example is the combination of fourth- and 
fifth-order Runge-Kutta schemes [198] . As a low-order example, we apply the 
outlined strategy to a second- and a third-order method. Let ^{t) be the exact 
solution of the differential equation and assume that an approximation to the 
solution at is known: ~ y"'. Take three RHS evaluations of the ODE: 



fi = f{t^,yn (90) 

/2 = /(r + At,y" + At/i) (91) 

f3 = f{t- + ^,y- + ^[h + M) (92) 

to produce a second-order (this is just the modified Euler method) 

a^^yn + At^A±Jl (93) 

and a third-order method 

a2 = |/'*+^[/i + /2 + 4/3]. (94) 

Thus, we have 

$(r + At) = oi + KAt^ + 0{At^) - 02 + 0{At^). (95) 
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Therefore, the leading error is 

E=\a2-ai\= KAt^ + O(Af^) (96) 

and the error grows at a rate 

^ ~ S ~ XAi^ (97) 

The step can either be accepted (e < e'^^'^), t"'^^ = t" + At, = 02, or 

rejected (e > e'^^'^) and re-taken. In both cases one can use the error measure 
to estimate the new time step. Since e cx At^, one finds 

At" , 

(98) 



gtried V /\y;tried 



and the suggestion for the new time step is 

(_acc \ 1/2 
^) ^^"^^^ (99) 



where s < 1 is a "safety factor" for a conservative choice of the next time step. 
For further reading on ODE-integration the excellent text books of [16,192,193,199] 
are recommended. 



2.9 "Best practice" suggestions 



We want to collect here a couple of suggestions born form practical experience 
that should help to carry out reliable simulations and to avoid numerical 
artifacts. This is somewhat "soft", yet hopefully still useful knowledge. 

• Neighbor numbers and resolution 

We recommend a large neighbor number, typically 100 or more. This some- 
what increases the smoothing lengths and thus deteriorates the integral 
approximation, Eq. (18). On the other hand, it can substantially reduce 
numerical noise, i.e. fluctuations around the true solution. To compensate, 
it is advisable to run simulations at the highest affordable particle number, 
numerical resolution does matter. Substructures can only be considered re- 
solved if they are substantially larger than the local smoothing lengths. If 
just an idea about the mass distribution is required, one may get away with 
a small particle number, say a few thousand particles, but reliable thermo- 
dynamic properties usually require much larger particle numbers. Contrary 
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to some prejudices, SPH can resolve shocks properly, but since shocks are 
spread across a few smoothing lengths an appropriate numerical resolution 
is required. 
Pzirticle masses 

If at all possible, equal-mass particles should be used. Light particles that 
interact with much more massive ones easily become "nervous" ("ping pong 
and cannon ball effect"), i.e. they can exhibit large fluctuations. In some 
problems different particle numbers are unavoidable, if, for example, a highly 
centrally condensed star is considered, the Courant time step criterion may 
restrict the central time steps to prohibitively small values. In such cases one 
should make sure that interacting particles only differ in masses by factors 
of very few. 
Initial conditions 

Accurate initial conditions are absolutely crucial for reliable simulations. 
The SPH particles should start out from their true numerical equilibrium 
positions. This can be obtained by applying an artificial damping term 
oc Va in the momentum equation that drags the particles to the desired 
initial conditions. Once equilibrium is found the particles should not move 
considerably during a few dynamical time scales after the damping has been 
switched off. 
Time integration 

Ideally, an integration scheme should monitor the achieved accuracy for the 
tried time step and reject it if the error becomes unacccptably large. Just 
estimating a time step that seems reasonable and than hope for the best, 
may produce noisy, or in the worst case, spurious results. The rejection of 
time steps in individual time step schemes may, however, pose bookkeeping 
challenges. A simple workaround has been proposed by [200]. For each and 
every simulation the conservation of energy, momentum and angular mo- 
mentum should be monitored. Reducing the time step size and increasing 
the force accuracy, say, if a tree is used for gravity, should improve the con- 
servation properties. A correct code should ensure conservation to better 
than 1% over several thousand time steps. 
Artificial dissipation 

Artificial dissipation should only be applied where really necessary. Schemes 
with time-dependent parameters are highly recommended. Modern versions 
of artificial dissipation terms such as Eqs. (67) and (68) should be used. 
Recent years have seen a large leap forward in eliminating unwanted effects 
from artificial dissipation, yet, there is certainly room for further improve- 
ment. 

Displaying results 

Generally, display continuous quantities and avoid particle plots to present 

physical results. The particles are just auxiliary constructs, "moving inter- 
polation points" that represent a continuum. Just showing particle positions 
projected on a plane may mimic a "resolution" that does not exist, a con- 
tour plot is a more honest presentation. To check that the numerics is well 



36 



behaved, it may, however, be useful to inspect the particle distribTition. A 
noisy particle distribution should always raise doubts about the quality of 
a simulation. 



Summary of "vanilla ice" SPH 

We summarize here the most basic form of the SPH equations. As the particle 
masses are kept fix, there is no need to solve the continuity equation. Densities 
can be obtained via summation from 

Pa^Y^^bWab. (100) 
b 

Alternatively, the continuity equation can be integrated, sec Eq. (3). The evo- 
lution equation for the specific internal energy can be written as 



du, 
~dt 



- = E^f + Vab ■ ^aWab (101) 



and the momentum equation as 

^ = -E^^f^ + 4+ VaWab. (102) 

b \Pa Pb J 

Particles are advanced in time according to this equation. Uab is the "standard" 
artificial viscosity term given in Eq. (59), more modern prescriptions are given 
in Eq. (67) and Eq. (68). Alternative forms of energy may be integrated forward 
in time, see, for example, Eq. (34). 



3 SPH from a variational principle 



In the previous section we had seen how the "vanilla ice" SPH equations can 
be derived directly from the Lagrangian form of the inviscid hydrodynamic 
equations. Although this form ensures conservation of energy, linear and angu- 
lar momentum and works well in practice, the corresponding symmetries were 
enforced somewhat ad hoc. SPH can also be derived using nothing more than 
a suitable fluid Lagrangian, the flrst law of thermodynamics and a prescription 
on how to obtain a density estimate via summation [165,201,146]. 
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3.1 The Lagrangian and the Euler-Lagrange equations 



The Lagrangian of a perfect fluid is given by [202] 

/o,2 \ 



P ( y - u{p, s) 1 dV, 



(103) 



where p is the mass density, v the fluid velocity, u the specific energy and s 
the specific entropy. In SPH-discretization it reads 



(104) 



The discretized equations for the fluid are then derived by applying the Euler- 
Lagrange equations 



d ( dL 



dL 



dt \dVa) dfa 



0, 



(105) 



where fa and Va refer to the position and velocity of particle a. The term 
in brackets yields the canonical particle momentum (like usual, we keep the 
particle mass fixed in time) 



dVa 



d_ 

dVn 



J2^b Y - 1i(P6, Sb) 



(106) 



and therefore, the first term in Eq. (105) yields the change of particle momen- 
tum rua^. 

The second term in the Lagrangian acts like a potential, and if self-gravity is 
considered, — nT-bUb has to be replaced by — ^b(wb + *^'fe); where is the 
gravitational potential [203]. In the following, we will only consider = 0. 
The second term in Eq. (105) becomes 



dL _ d 
dfa dfa 



u{Pb, Sb) 



dub 
dpb 



dpb 
dfa 



(107) 



We can make use of the first law of thermodynamics, Eq. (7), to find 
dva ^ Pb dpb 



m. 



- ^-^rub- 

b 



dt V Pb dVa 



(108) 
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3.2 The density, its derivatives and "grad-h" -terms 



So far, we had ignored all extra terms that result from variable smoothing 
lengths [204,146,147], this is what we address now. Prom now on, we use the 
smoothing length of the particle itself in the density summation, 

Pa = J2^bWirab,ha). (109) 
b 



and use Eq. (51) to adapt the smoothing length. Thus, pa depends on ha and 
vice versa, which requires an iteration to reach consistency. 
If we take the changes of h into account, the Lagrangian time derivative of 
the density is given by 



dpa d rrr f , ^\ ^ i dWab{ha) drab , dWab{ha) dha 



^ dWabiha) . ^ , V- dWab{ha) dha dpa 

= E rrib^—eab ■ t'a^ + E "^'> " ^ ^ 

^ ^ T.r ^ .dhadpa ^ dWab{ha) , . 

= E m.Vab ■ VaWabiha) + E^^^^^ (HO) 

where we have used Eqs. (24) and (26). If we collect the dpa/dt-terms and 
introduce 

O - A 9ha^ dWab{ha) \ 



the time derivative of the density reads 

^ = ^T^rUbVab ■ VaWabiha). (112) 
dt i^aT 



This is the generahzation of the SPH expression for the density change that 
follows from Eq. (31). In a similar way, the spatial derivatives can be calculated 



dp 



b )v7 ui \^ dWbkih) dh dpb 



n 



b k 



(113) 
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To summarize: if derivatives of the smoothing length h are accounted for, the 
"standard" SPH expressions for the density derivatives have to be corrected 
by factors l/fl, see Eq. (Ill), (112) and (113). 



3.3 The SPH equations with "grad-h" -terms 



We keep again the masses fixed in time, so there is no need to solve the 
continuity equation. By inserting Eq. (112) into Eq. (30) we find the energy 
equation 

1 P 

For the momentum equation we need to insert the density gradient, Eq. (113), 
into Eq. (108) 



T^J^ = - y^rni^^WaPh = ( ^y^rrikW aWbkih,,)] . (115) 

Using Eq. (25), the above equation becomes 



ma^^-J2^^b^7^1Z'^k'^bWkb{hb) {Sba " Ska) 



dt ^ pI Clb 

= -ma^-^^rrikV aWka{ha) + ^rnb-^-^maVbWab{hb) 

Pl^ak b Pb 

PI PI 

Pa^^'a b b Pb^^b 

^-ma^rUb l^VaWabiha) + aWab{hb)] . (116) 

Thus the final momentum equation reads 

f = -E™. + ^VM)) . (117) 



Note that now all arbitrariness (such as the particular form of Eq. (28)) has 
been eliminated from the derivation. The "grad-h" terms increase the accuracy 
of SPH and the conservation properties in the presence of varying smoothing 
lengths. How important they are in practice depends on both the problem and 
the numerical resolution [146,112]. 
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SPH equations from a variational principle 

With 

(TTl \ 
— (118) 
Pa J 

the energy equation reads 
dv 1 P 

= 5: • VaWa^iK) (119) 

^^a Pa b 

and the momentum equation is 

^ = - E (TP^'^aWabiha) + aWabih)] , (120) 

dt Y V^aPa ^bPb J 

where 

(121) 

are the so-caUed "grad-h-terms" . Additional artificial dissipation terms as dis- 
cussed in Sec. 2 need to be applied in order to resolve shocks. 



4 Relativistic SPH 



Several relativistic versions of SPH exist. Early formulations due to Mann 
[205,206] were able to capture the overall solutions, but they were -by today's 
standards- not particularly accurate. 

Laguna et al. [207] developed a 3D, general-relativistic SPH code. Due to 
their choice of variables their continuity equation contains a gravitational 
source term. This term requires to introduce SPH kernel functions for curved 
space-times which are in general not isotropic and invariant under transla- 
tions. Moreover, they did not use a conservative form of the equations and 
their approach requires additional time derivatives of Lorentz factors which 
they approximated by first order finite differences. Similar to the approach of 
Mann, only mildly relativistic shocks could be treated. The approach of [207] 
is also the basis of the more recent code described in [208]. 
Kheyfets et al. [209] suggested an alternative approach in which they define 
the spatial kernel function in a local rest frame and assume that space is ap- 
proximately flat in this frame. The interaction between particles requires the 
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choice of a particular frame whose calculation requires the solution of foTir 
coupled algebraic equations. While this approach is completely covariant, it 
requires a large computational effort. 

Like in grid-based methods, a conservative formulation is crucial to handle 
strong shocks. As we had seen in Sec. 2. 7, concepts of Riemann solvers were 
used as a guide to improve SPH's artificial viscosity by treating the interac- 
tion between two particles as a local Riemann problem [183]. This approach 
replaces wave propagation speeds by an appropriately chosen signal velocity. 
With this approach (together with evolving the total energy rather than the 
thermal energy) Chow and Monaghan [210] obtained accurate results even for 
ultra-relativistic flows. 

More recently, Siegler and Riffert [211] and Siegler [212] have suggested a set 
of equations for both the special- and general-relativistic (fixed background 
metric) case. It is based on a formulation of the Lagrangian hydrodynamics 
equations in conservative form. By a particular choice of dynamical variables 
they avoid many complications that have plagued earlier rclativistic SPH for- 
mulations. In particular, by choosing the relativistic rest mass density for the 
SPH formalism the continuity equation has the same form as in non-relativistic 
SPH, time derivatives of Lorentz factors do not appear and flat space kernels 
can be used. This comes at the price of inverting the dynamically evolved vari- 
ables into the physical variables. Siegler [212] was able to accurately simulate 
some test cases with Lorentz factors up to 7 = 1000. 

A more elegant approach is based on the use of the discretized Lagrangian 
of a perfect fiuid [201]. Guided by the canonical momentum and energy, suit- 
able numerical variables can be chosen which lead to evolution equations that 
are similar to the Newtonian ones. Due to the consistent derivation from a 
Lagrangian the conservation of mass, energy, linear and angular momentum 
arc guaranteed. This approach can be applied both to the special- and the 
general-relativistic case. 

In recent years, SPH has also been applied to study flows in time- varying 
space-times. The first approaches used Post-Newtonian approximations [213- 
216], more recently the conformal flatness approximation [217,218] has been 
implemented [219-221]. 

In the following subsections we derive the special- and general-relativistic SPH 
equations from a variational principle similar to [201]. 

4-1 Special-relativistic SPH 

Here, we generahze the approach of Monaghan and Price [201] to include the 

special-relativistic "grad-h" terms. We use the conventions that the metric 
tensor has the signature (-,+,+,+), Latin indices run from 1 to 3, Greek ones 
from to 3 with the zero component being time. In addition, the usual Einstein 
sum convention and c — G — 1 are used and the flat spacetime metric is 
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denoted by rj/j,,,. 

With these conventions, the hne element reads 

^ rj^^dx'^dx'' (122) 



and the proper time is given by dr"^ — —ds^. The 4- velocity is defined as 

C/'^ = ^ (123) 
dr 

and, due to the above two relations, it is normalized to —1: 
dr^ dri^ r^s^ 

The velocity components can be written as 

{U^) = {l,iv), (125) 



where 7 is the Lorentz factor 

with V being the 3-velocity. Due to Eq. (125) the 3-velocity can be expressed 
as 

The 4-momentum of a particle with rest mass m is given by 

(p^) ^ {mU^) = {-fm,-fmv) = {E,p), (128) 

where E is the particle energy and p its relativistic momentum. The last 
expression is also correct for particles with vanishing rest mass. 

4.1.1 The Lagrangian 

Consistent with our overall strategy, we derive the SPH equations from a 
Lagrangian to ensure exact conservation. We start from the Lagrangian of a 
perfect fluid [222] 

i^pf,sr = - / Ti^^U^U, dV. (129) 
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The energy momentum tensor of a perfect fluid 

The energy-momentum tensor, T^*^, is a key element of relativistic hydrody- 
namics as it contains the sources of the relativistic gravitational field. Its com- 
ponents T'^" are defined as "flux of //-component of the 4-momentum in in- 
direction" . Thus, the T°°-component is the energy density, the T°-^ component 
is the energy flux, the component is the momentum density and the T*-^ is 
the momentum flux. For a perfect fluid with zero viscosity and heat conduction 
it is given by 

T^"" = (P + e)C/'^C/^ + Pt)^'', (130) 

where P is the fluid pressure, the 4- velocity and e the energy density in the 
local rest frame. The conservation of energy and momentum can be expressed 
as 

T^^^/ = 0, (131) 

where the comma denotes the partial derivative (in general relativity it has 
to be replaced by the covariant derivative, usually denoted by a semicolon). 
11 = expresses energy and = i momentum conservation. 
In the Newtonian limit (with c = 1) the following relations hold C/° = 7 ^ 1, 
[/* ~ V*, e ~ p and P/e <^ 1. Thus the components read 

T°° = (P + e)t/°f/° + P^p (132) 
T°J' = (P + e)U°W ^ pv^ (133) 
T'^ = (P + e)U'W + P6'^ fti pv'v^ + PS'^. (134) 

Therefore, T^^,i^ = reduces to the usual continuity equation 

= T'',o +T'^,j = 1^ + ^ = (135) 



and T^^,n = becomes 

_ ^ 

dt ' dx^ ' dx^ 



r^,^T%+r-^%;i + '^ + |5 = o, (136) 



which is the Euler equation. 

This Lagrangian is, apart from the fiat space volume element, the same as the 
general-relativistic one, see Sec. 4.2. The energy density, , that is needed in 



^ This quantity is usually denoted by p in the Uterature. To avoid confusion with 
the non-relativistic mass density that we have used previously, we choose a different 
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the energy-momentum tensor, see Eq. (130), has contributions from the rest 
mass and the thermal energy: 

e = Brest + etherm = PrestC^ + ^^Prest = nmoC^{l + u/c^), (137) 

where mo is the baryon mass and u is again the specific thermal energy and 
n is the baryon number density, the latter two are measured in the local rest 
frame. Here, we have kept the c^s for clarity. From now on we will measure 
all our energies in units of moc^ (and use c = 1). With these conventions the 
energy momentum tensor reads 

T^'^ = (n[l + u{n, s)] + P)U^U'' + (138) 



Here s is the specific entropy per baryon in the local fiuid rest frame and P the 
fluid pressure. For a practical computation we flx a particular frame ("com- 
puting frame") and transform the quantities into this frame where necessary. 
Using Eq. (124), we can considerably simplify the above Lagrangian. The 
quantity under the integral becomes 

T^'U^U, = {(n[l + u\ + P)UW} U^U, + Pv'^'U^U, = n{l + u) (139) 

and the Lagrangian simplifies to 

V,sr ^-Jn{l + u) dV. (140) 



The baryon number density in the rest frame of a fluid element, n, is related 
to the number density in our computing frame, A^, by a Lorentz factor, 7. To 
see this, assume a volume in the local rest frame of a fluid clement, 1^^, that 
contains k baryons. In general, this fluid element will move with some velocity 
with respect to our computing frame. If we assume that the motion is in, 
say, X— direction, the x— component of the volume element in the computing 
frame, Kf, appears Lorentz-contracted by a factor 7: V^f = Vm/i- Since this 
volume element contains the same number of baryons, k, the density in the 
computing frame is 

(Ml) 



As before, see Eq. (104), we want to discretize the Lagrangian to derive the 
SPH equations. We subdivide space in volumes AV^, so that a volume labeled 

symbol here. 

Note that this quantity depends on the relative number of neutrons and protons, 
i.e. on the nuclear composition. 
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by index b contains ui, ~ AVbNi) baryons, where is the baryon density at 
particle b (all quantities measured in the computing frame). Or, the other way 
around, if a fluid parcel labeled by b contains Uf, baryons, it has a volume 

AH = ^. (142) 



Similar to Eq.(12), we can use this volume element in the SPH-approximation 
of a quantity / : 

m-T.fb^W{\r-n\,h). (143) 



b 



With this prescription the baryon number density in the computing frame at 
position of particle a reads 

Na = N{ra) = Y.l^bW{\fa - n\,ha), (144) 
b 



where we have used ha to conform with the non-relativistic version, Eq. (109). 
The latter can be recovered by the following replacements: — ^ z/5 and 
P6 Nil. If we keep the baryon number per SPH particle, Uh, fixed, the 
total baryon number is conserved and there is no need to evolve a continuity 
equation. Analogously to Eq. (31) one could also discretize the continuity 
equation and calculate the density by integration. We adapt the smoothing 
length similar to the Newtonian case, see Eq. (51), 

K = r,{£ . (145) 

which, again, requires an iteration for consistent values of Na and ha- We can 
now discretize the Lagrangian of Eq. (140) 




(146) 



or, by using Eq. (141), 

i^SPH,sr = - E + Sb)] =-J2^b - ^b [1 + (147) 

b h 



For ease of notation we have again dropped the distinction between the function 
to be interpolated and the interpolant. 
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4-1.2 The momentum equation 

The Euler-Lagrange equations, Eq. (105), determine the evolution of special- 
relativistic momentum: 

Pa=^^ ^-l^'^bjr^l • (148) 



The specific energy u depends, via the density, on the velocity: 

dub f dub \ dub Pb drib 
dva [dubldva nldva 



(149) 



where we have used Eq. (8). Via Eq. (141) we have 

^ = Nb^ (-] = Nb^{l - vlf' = -NbibVb5a,. (150) 

dVa dVa V76/ 

Thus, the canonical momentum is 

Ub du 



IbdVa 



T^b ( Pb 



= ^b{VblbKb)[^ + ^^b] + H — ^ (Nb-fbVbSab) 
b b 76 

= J^alaVa f 1 + + — ) , (151) 

where we have again used Eq. (141), the term in brackets after the last equal 
sign is the enthalpy per baryon. Obviously, in relativity the pressure con- 
tributes to the momentum density. We define the canonical momentum per 
baryon as 

Sa = ^aVaU + Ua+—), (152) 



whose time evolution is governed by dLsp}i,sr/dfa, see Eq. (105), and which 
requires the gradient of the number density. The latter is obtained similar to 
the Newtonian case, see Eq. (113), 



dNb_y^^ j dWjr-bk, hb) dnk ^ dWjrbk, K) dhb dNb ' 
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i E ^kyaW,k{K){5u - ha) (153) 



with 



dhb ^ dW{rbk,hb) 



Thus 



^-^SPH,sr _ sr^ l^b dUb _ UbPbdUb _ Vb Pb dNb 

dVa b 76 dfa Y 76 b 76 ^6 ^^a 

^ ~;^VaW^6ik(/l6)(56a " (^fea) 

% ^hrib 



= -i^a E ^b I V„W^„6(/l„) + ^ V„W^„6(/l6) \ , (155) 



where we have used Eqs. (8), (141), (144), (153), and (26). With Eq. (155) 
and Eq. (151), the special-relativistic SPH momentum equation becomes 

^ = - E ^ J T^.'^aWabiK) + ^,VaWab{hb) ] • (156) 
b l^a^a ^bN^ J 



This form can be obtained from the non-relativistic case, Eq. (117), by the 

— * 

following replacements: Va Sa, ruk — > ffe, ^k and pk — > Nk- 

Since the Lagrangian is invariant under infinitesimal translation the total 

canonical momentum 

P = y^Pb = E ^bSb = E ^blbVb (l+Ub+ — ) (157) 



is conserved. Similarly, due to invariance under rotation, the angular momen- 
tum 

L = E n> X Pb = E X Sb = y] i^blb (l + Ub+ — ) fb X Vb (158) 
K u V V UbJ 



is conserved by construction. 
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4-1.3 The energy equation 



To choose a suitable numerical energy variable we start from the conserved 
canonical energy (no explicit time dependence in the Lagrangian) 

E = 2_^ ■ Va - LsPH,sr- (159) 



With Eqs. (147) and Eq. (151) the energy reads 



E^^Ua [va ■ Sa+ llt^ j = ^ 
a \ la / a 



(160) 



where, we have defined 



ea = Va- Sa+^—^^. (161) 
7a 



For later use, we note that by using vl — 1 — 1/7^ and Eq. (141) we can 
express as 

To find the evolution equation for we need the time derivative of the second 

term in Eq. (161). Ideally, it should not contain a time derivative of the Lorentz 
factor since such terms are known to destabilize numerical schemes [223] . Using 

^=7„^a-^ and -^ = Jt[-)=-^-Na7aVa-^{lQ3) 

we find 



a 



d fl+Ua\ if dUa , .C?7, 

1 Pa dUa 1 + Ua 3 ^ dVa 
-TaVa ■ 



7a nl dt 72 " dt 

PadNa A , Pa\ ^ d.V„ 



1 + Ma H 7at'a 



m dt V rinJ ' dt 



m dt dt ' ^^^^^ 



a 
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where we used Eqs. (8), (141) and (163). With equation (164) at hand, we can 
take the derivative of Eq. (161) 



— - — l-y -S + lltJ^l 
dt dt\"'"' la \ 

—* 

~ dt dt ^7V2 dt dt 

— * 

^ dSa , Pa dNa 

i.e. apart from the relativistic momentum equation (156), we need the La- 
grangian time derivative of the number density in the computing frame: 



dNa_y. f dWabihg) drab dWab dhg dNg 
dt drab dt dha dNa dt 



^-^Y. ^bVab ■ ^aWab{ha), (166) 
"a b 

where we have used Eq. (154). If we insert Eqs. (156) and (166) into Eq. (165) 
the energy evolution equation becomes 

f = - 1: (H ■ V.H^.O.) + g ■ VM)) ■ (167) 



This relativistic energy equation looks similar to the non-relativistic equation 
for the thermokinetic energy, see Eq. (34). 



4-i-4 Recovery of the primitive variables 

Due to our choice of the numerical variables neither the momentum nor the 
energy equation contain time derivatives of Lorentz factors. Such terms re- 
stricted earlier versions of relativistic SPH to moderate Lorentz factors only. 
But this comes at the price that the physical variables need to be recovered 
at every time step from the numerical ones, a price that relativistic Eulerian 
codes also have to pay. 

At the end of a time step 7/f , u and n need to be calculated from the updated 

— * 

numerical quantities N, S and e. One can express all variables in the equation 
of state 

P„ = (r - l)naUa (168) 
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Fig. 8. Mildly relativistic (7max ~ 1-4) shock tube test: comparison between the 
Newtonian (dashed) and the special-relativistic results (solid, red). Upper left: ve- 
locity in units of c, upper right: thermal energy, lower left: pressure, lower right: 
(computing frame) number density N and Newtonian mass density p. 

as a function of the updated numerical variables and the pressure Pa- The 
resulting equation is solved numerically for Pa and once it is found, the other 
physical variables can be recovered [224,183]. From Eqs. (152) and (162) one 
finds 



ea + Pa/Na 



(169) 



and thus the Lorentz factor is 



7a 



l-S^Jiia + Pa/Naf 



(170) 
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Prom Eq. (169) we have 




Ian, 



a 



'a 



] 



Va^ JaVa ( 1 + l^a + 



Pa 



n, 



'a 



) 



(171) 



which can be solved for 



7a laNc 



r(l-7a')-l- 

a 



(172) 



With aid of Eqs. (141) and (172), Eq. (168) can be solved, e.g. via a Newton- 
Raphson scheme, for the new pressure Pa- Once Pa is known, the Lorentz 
factor can be calculated from Eq. (170), the specific energy from Eq. (172) 
and the velocity from Eq. (169). 



4.1.5 Numerical tests 

To explore the performance of this SPH formulation we show a numerical test 
that can be considered a special-relativistic generahzation of Sod's shock tube 
test, sec Sec. 2.7. This test has become a widespread benchmark for rclativistic 
hydrodynamics codes, see for example [224,210,211,225,226]. The test is per- 
formed with a polytropic equation of state with exponent P = 5/3, vanishing 
initial velocities everywhere, the left state has a pressure Pl = 40/3 and a den- 
sity Nl — 10, while the right state is prepared with Pr — 10~^ and Nr — 1. 
Although the resulting velocities are only mildly relativstic (7max ~ 1-4), the 
deviations from the purely Newtonian result are already substantial. This is 
demonstrated in Pig. 8, where we compare the exact solutions for these initial 
conditions for the Newtonian (dashed) and the special-relativistic case (solid). 
Pig. 9 shows the SPH result (black circles, about 3000 particles) at t= 0.35 
together with the exact solution (red, solid line). The general agreement is 
excellent, noticeable deviations are only observed in the form of a slightly 
smeared out contact discontinuity at x ^ 0.25. A striking difference to earlier 
SPH formulations [207,211] is the absence of any spike in u and P at the con- 
tact discontinuity. 

This SPH formulation has been further explored in a large set of special- 
relativistic benchmark tests [227]. As expected, it performs very well in pure 
advection problems. Maybe more remarkable, it also shows convincing results 
in extremely strong, rclativistic shocks. Por example, it yields accurate solu- 
tions in a wall shock test with a Lorentz factor of 7 = 50 000. The special- 
relativistic "grad-h" -terms generally improve the accuracy, but (at least in the 
performed set of tests) only to a moderate extent. Por more details we refer 
to [227]. 
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Summciry of the special-relativistic SPH equations 

The momentum equation derived from the special-relativistic Lagrangian reads 

dSn. 



dt 
where 



JC ( P P ^ 



5„ = Tati'„(l + «„+— ) (174) 



and 

The (canonical) energy variable 

ia = V, ■ Sa + (176) 

7a 

is evolved according to 

^ = - V Z/, f ■ VaWaM + • VaWabih)] . (177) 

The computing frame number density can be found either by summation, 

Na^^l^bWaM, (178) 
6 

or, alternatively, by integration of 

^ = 7^ E '^bVab ■ VaWM- (179) 

dt ^ 



4-2 General-relativistic SPH on a fixed background metric 



4.2.1 The strategy 

We now generalize the previous approach by assuming that a fixed metric 
is given as a function of the coordinates and that corrections to the metric 
induced by the fluid are negligible. Such a situation occurs, for example, during 
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Fig. 9. SPH solution of the relativistic shock tube test of [224] (shown are 3000 
particles; upper left: velocity in units of c, upper right: thermal energy, lower left: 
pressure, lower right: (computing frame) number density N). The black circles show 
the SPH solution, the red line marks the exact solution [226]. 

the tidal disruption of a star by a supermassive black hole in the center of a 
galaxy. Again, we start from the discretized Lagrangian of an ideal fluid and 
use the canonical momentum/energy as a guide for the choice of the numerical 
variables. As a matter of course, this new set of equations should reduce in 
the flat-space limit to the special-relativistic equations, see Eq. (173) - (178), 
which, in their low-velocity limit, should be equivalent to the non-relativistic, 
standard SPH equations. 

4-2.2 The Lagrangian 

We start from the general-relativistic Lagrangian of a perfect fluid [222] 

i^pf,GR = - / T^^'U.U, dV, (180) 

where dV is the proper volume element and g = det{g^^) is the deter- 

minant of the metric tensor, ^f^y. As before, the energy momentum tensor is 
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defined by 

T^'^ = {n[l + u{n, s)] + PjU^'W + Pg^"", 



(181) 



where the fluid quantities are measured in their local rest frame and the four- 
velocity is defined as in Eq. (125). Obviously, 

=(-*..'V)S (182) 



where we used — dx'^ /dt. Comparing with the special-relativistic relation 
dr/dt — 1/7, one sees that the quantity 

e = (-5^.^;'^^;^)-5 (183) 



takes over the role of and reduces in the flat spacetime hmit to the special- 
relativistic Lorentz factor 7 = — v"^. Similar to Eq. (127) we can also 
write 

„ dx^" dxf'dr W , 



where we have used dt/dr = = O, see Eq. (182). Using the normalization 
condition of the four- velocity, Eq. (124), the Lagrangian reduces to 



V,GR = - / [{n{l + u) + P}U>'U''U^U, + Pg>"'U^U,] y/^ dV 

= - y n[l + u{n, s)] dV, (185) 



which is, of course, up to the volume element, the same result as in special 
relativity, Eq.(140). 

To discrctize the Lagrangian, we flrst have to flnd a suitable density variable 
so that we can apply the usual SPH discretization prescription. The general- 
relativistic baryon number conservation equation reads 

^ ^- {V=^U^n) = 0, (186) 



—g dx^ 

see e.g. [228], and this suggests to use the modified number density 

N* = C/° n = e n, (187) 
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since it can be written with the help of Eq. (184) as 

ON* d(N*v') , , 

Eq. (187) generahzes the carher relation between the computing frame and 
the local rest frame, Eq. (141). We use N* in the SPH interpolations 

/(0=EA^^(I^"'-^1>I)> (189) 



where 

= iV*AH = nfoGfev^AH (190) 

is the number of baryons contained in the volume y/—gbAVb. Applying Eq. (189) 
to N* one obtains a density estimate by summation 

K = T.^bWab{ha), (191) 
b 

similar to the earher results, Eqs. (109) and (144). We can now use the relation 
between N* and n, Eq. (187), to write the Lagrangian as 



Lpi,GR = -] ^=jQ [1 + u{n, s)]y^dV = -J—[l + u{n, s)]dV, (192) 

i.e. the metric tensor contribution ^/—g drops out, and using AVb — Ub/N^, 
the SPH form of the Lagrangian reads 

-f'SPH.GR = - 7^ [1 + '^i^b, Sb)] ■ (193) 
h '^b 



In the flat spacetime limit ©5 — jj, and we recover the discretized, special- 
relativistic Lagrangian, Eq. (147). 



4.2.3 The momentum equation 

To identify a suitable numerical variable, we use like before, see Eq. (152), the 
canonical momentum 



dL 
dv 
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We thus need 

4 ^ 4^-^-^'^^)^'' = ^9^,vna So, (195) 



and the derivative dub/dvl, where, hke in the special-relativistic case, the 
velocity dependence comes from the Lorentz-factor-hke quantity © relating 
the densities in the local frame and the computing frame 



dub dub dnjj Pb d 



Pb d 



dvi dubdvi nldvl\^/^^Qb, 



nt 



-9b 



dvi 



1 

©6 



Pb m 



nt 



-9b 



-Qbi9inV^)a Sab- 



(196) 



Here we have used Eqs. (8), (187) and (195). By inserting Eqs. (195) and (196) 
into Eq. (194) we find 



dvi 



^6[1 + Ub] {-Qb {9iX)a Sab) - E ^ 
b b 

I P N* \ 

Uaea{l+Ua+^ ^ ) (5iX)a 



■„e„ (^1 + + (yiX)a' 



Pb N* 

^ Qfc {giluV'^)a Sa, 



H \/-9b 



V, 



(1^ 



where Eq. (187) has been used. This is the i-component of the canonical 
momentum and, like before, we use the canonical momentum per baryon. 



^ = eJ 1 + Ii„ + ^ {9,^V^)a. 

Vadv\ \ UaJ 



(198) 



as numerical variable. Again, the term in the first bracket is the specific en- 
thalpy. In the fiat-space limit, ©„ ja, g^u ^ Vnu and gif.v'' ^ Vi = v\ 
this expression reduces exactly to the special-relativistic case, Eq. (174), as it 
should. 

The evolution of Si^a is determined by the Euler-Lagrange equations, so we 
need 



The first derivative can be written as 
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d 1 



d 



e,fdg^\ ^ ^^^^ 



2 V dx 



a / b 



(200) 



and due to Eq. (184) 



02 ' 



(201) 



so it becomes 



I 26 )i 



.9x1 J, 



(202) 



The derivative of the internal energy is 



duh duh driK Ph d 



dxi 



drib dx'a 
Pb 



nl dxi 
(dNf 



^b ©6 



-g. Qb V dxi 



+ 



PbNl 



d 1 



nt 



b^b \(^K^b, 



PbN; 



d 



nlQb \dxi./=g^^ 



(203) 



where we have used the first law of thermodynamics, Eq. (8) . By using Eq. (25) , 
the derivative of the number density becomes 



dN* _ d 
dxi 



a \ k J k "^b 



(204) 



The last remaining derivative is 



d 1 



g"' dg^f 



9x1 \F9b 



2^=^ dx^ 



^ab, 



a / b 



(205) 



see the box on the derivatives of the metric tensor determinant. 
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Derivatives of the metric tensor determinant g 

We will collect here some formulae related to the quantity g that are needed 
in several places. The derivative of the metric tensor determinant is generally 
given by (see e.g. [228]) 



dg ^ a^ dggf) 



Therefore 



2 ^ dxi' 



_d_ 
dx' 



9, 



jr_dg^ 
2^=^ dx^ 



(206) 



(207) 



and 



(^) = ^ (^) 



dt 



dx^ 



^ ap9gal 

2 ^ 



-9 „al} dgal} 



dt 



(208) 



Of course, if the metric tensor at the position of particle b is to be differentiated 
with respect to a property of particle a, the appropriate Kronecker deltas, 5ab, 
have to be applied. 

With the help of Eqs. (202), (203), (204) and (207) Eq. (199) becomes 



dL 
dxi 



26 



dx 



Jab 



a / b 



PbNl 



2^^k^^{dba-dka) 



+ 



+ 



,2^ 

PbNl 
nlOb 



dxl 

V 20 )b\dx\^ 
9^" 



^ab 



-g dxi 



^ab 



(209) 



The terms that involve the derivatives of the kernel represent the hydrody- 
namic part of the equations, the others the action of gravity. After eliminating 

one sum in the hydrodynamic terms via the Kronecker symbol, relabeling the 
summation index from k to b and using kernel property Eq. (26) they read 



'dL\ 



-T^a 1^ T^b ( 2 n2 + 



ab 



nt 



(210) 



and on using Eq. (187) 
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dL 



V a/ h 



]\T*2 
a 



+ 



dxi ' 



(211) 



The terms that involve derivatives of the metric ("gravity terms") can be 
written as 



dxi 



2e„ 



P N* \ P 



ni 



(212) 



If we apply once more Eq. (187) and remember the form of the energy- 
momentum tensor of the fluid, Eq. (181), the above term simply becomes 



\dxi 2N* V dx' ' 

Putting all together, the spatial derivative of the Lagrangian reads 



(213) 



dL 

dxi 



I b 



*2 



+ 



9bPb\ dW, 



ab 



N, 



*2 



dxi 



+ 



9a 



2N* 



dg, 



dx' 



.(214) 



Note that the first, hydrodynamical term has the symmetry in the particle 
indices, a and b, that we know from previous forms of SPH equations, while 
the gravity term is exclusively evaluated at the position of the particle under 
consideration, a. 

Inserting Eqs. (197), (198) and (214) into the Euler-Lagrange equations yields 
the final, general-relativistic momentum equation 



dt 



b 



-9aPa 



]\T*2 
a 



+ 



-g,Pb\ dW, 



ab 



dxi 



+ 



-9a (rj.^,ud9t.u 



2N* 



dx^ 



.(215) 



4.2.4 The energy equation 

We use once more guidance from the canonical energy to choose the numerical 
energy variable: 



a '^^a a "a 



(216) 
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Here we have introduced the energy per baryon 



■Qail + ua + ^) {9^Xv% + ^ = ^.X + (217) 

very similar to the special-relativistic energy variable, Eq. (161). Its temporal 
change is given by 



dia dSi,a i d f l + Ua \ 

The first time derivative is known from the momentum equation, the second 
can be hoped to cancel out with a term resulting from the last term, as in the 
special-relativistic case, see Eq. (164). For this last term we need dua/dt and 
dQa/dt. To obtain the time derivative of Ua we use the first law of thermody- 
namics, Eq. (8), the number density relation between the frames, Eq. (187), 
and d{y^)/dt, Eq. (208). This delivers 



dua Padf N* \ P^^^Q^dN: ^ P,N: d f I \ 



dt nldt\^/=gej^ Nf dt nl dt \y/=gQ ) ^ 

The last derivative is given by 

d_ (^_\ ^ _ (^_ .,d_g^\ _ (^_dQ\ 
dty^Qj^ [2^Q^ dt y^gQ^dtjJ ^ ^ 

thus the change in Ua becomes 



dUa _Pa^/^a^CidN* Pa^^a^dQc^f} Pa dQa 



aji "'^ciP _ - a a (221) 

dt Nf dt 2n„ dt n„e„ dt ' ^ ' 

It remains to calculate the derivative of the generalized Lorentz-factor 

On using Eqs. (221) and (222) we find 



d n 
dt \ 



9a^a dNl 




Pg a0 dgafS _ Pg dQ 

2na dt UaOa dt 



1 + Ua dQ 

02 lit 
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-9n dN* 

J a a 



]\T*2 
a 



dt 



-ea[l + Ua + 



q"'" H 



l + Ua + 











naJ 


Pa 




ria- 





lit J 

dgap 
dt 



(.223) 



If we replace the velocities via Eq. (184) by four- velocities and use the defini- 
tion of the energy momentum tensor, Eq. (181), the equation simplifies to 



d^ / I + 
dt \ Ba 



df_' 
dt 



2N* dt 



gT<-Pdgo^p 



PgyTg, dN: dvl_ 

N*^ dt ^ dt \ 2N* dt 



(224) 



Here we have split the contraction in the term containing dv^ /dt into a spatial 
and a temporal part and we made use of = 1, see Eq. (184). This equation 
is very similar the special- relativistic one, Eq. (164), but it also contains a 
gravity contribution. If we insert Eq. (224) into the energy evolution equation 
(218), use the SPH prescription for dN* /dt and collect terms into "hydro" 
(kernels) and "gravity" (metric derivatives) we find 



dt 



b 
b 



-dgPaVb 

a 

^9aPaVb 

]Vf*2 
a 



+ 



+ 



2N*~ " 



This is the final, general-relativistic SPH energy equation. The equations 
(191), (215), and (225) together with an equation of state form the set of 
general-relativistic SPH equations in a fixed background metric. 
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Summciry of the general-relativistic SPH 
equations on a fixed background metric 

Ignoring derivatives from the smoothing lengths, the momentum equation 
reads 



where 



S,,a = Qa(l + Ua + — ) {g^X)a (227) 



is the canonical momentum per baryon and 

Qa = {-g^.v^vX'' (228) 
the generalized Lorentz factor. The energy equation reads 

f = - E + ^«-») • V.l^. - ^ (r-^) .(229) 



where 



ia = S,,aVl + (230) 

is the canonical energy per nucleon. The number density can again be calcu- 
lated via summation, 

K = ll^bWab{ha). (231) 



5 Summciry 



Since the Smooth Particle Hydrodynamics method was suggested in the late 
seventies it has undergone a long sequence of technical improvements. In par- 
allel, the method has been applied to a large variety of problems both inside 
and outside astrophysics, and consequently a slew of different physical pro- 
cesses has been included in SPH-based simulations. 

In this review we have only very briefly touched upon this latter development, 
we merely provided pointers to the current literature. Instead we have focused 
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on what we hope is a pedagogical introduction leading to an in-depth under- 
standing of how the method works. All essential equations were derived, and in 
particular, a commonly used ( "vanilla ice" ) set of SPH equations that directly 
discretizes the Lagrangian equations of an ideal fluid. We have also reviewed 
concepts such as adaptive numerical resolution, the reasoning behind artificial 
viscosity and basic aspects of the ODE integration in the SPH context. 
The "vanilla ice" equation set works well in practice and possesses the property 
of "hard-wired" conservation of mass, energy, linear and angular momentum. 
The symmetrization in the particle indices is, however, somewhat arbitrary 
and was enforced ad hoc. In Sec. 3, a modern form of the SPH equations was 
derived that improves on this weakness: it uses nothing more than the La- 
grangian of an ideal fluid, the first law of thermodynamics and a prescription 
how to obtain the density by a summation over particles. After the Lagrangian 
has been written in its SPH-discrete form, there is no more arbitrariness left 
in the rest of the derivation: the equations follow stringently from the Euler- 
Lagrange equations. In addition to curing the aesthetical shortcoming of the 
previous approach, the latter also naturally leads to corrective terms due to the 
derivatives of the smoothing lengths, in SPH usually called "grad-h-terms" . 
The elegant variational concept naturally carries over to both the special- and 
general-relativistic (fixed-metric) case. Here, we derived for the first time the 
special-relativistic SPH equations that include "grad-h-terms" . The numerical 
variables are chosen as canonical momentum and energy per baryon. While 
this has obvious numerical advantages, it comes at the price of recovering the 
physical variables from the numerical ones at each time step, a tribute that 
is also payed in modern grid-based relativistic approaches, see e.g. [226] . The 
performance of this new equation set has been briefiy illustrated at the exam- 
ple of a relativistic shock tube test, for a more exhaustive set of tests we refer 
to [227]. 

We conclude this review with a derivation of the general-relativistic case in 
which the space-time metric can be considered unperturbed by the self-gravity 
of the fiuid, a case that is realized, for example, in the tidal disruption of a 
star by a supermassive black hole. The resulting equations are a natural gen- 
eralization of the special-relativistic case and do not introduce much of an 
additional numerical burden. The resulting extra terms can be analytically 
calculated from the metric that is assumed to be known. 
In the future, more physical processes and their corresponding numerical 
schemes will certainly find their ways into SPH. On the numerical hydrody- 
namics side the recent past has seen a couple suggestions on how to improve 
on known weaknesses of SPH. One of them was SPH's unsatisfactory per- 
formance on fiuid instabilities across discontinuities with large density jumps 
[185]. A modification of artificial dissipation terms [154] or, alternatively, a 
more general formulation of the SPH equations [162] have been suggested 
as cures for this problem. In terms of artificial dissipation, a tensor artificial 
viscosity approach such as the one suggested in [175] should definitely be ex- 
plored further for its potential and usability in astrophysics. While the special- 
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relativistic SPH equations as derived in Sec. 4 have been carefully tested, the 
general-relativistic equation set has so far not been implemented and applied 
to astrophysical problems. This task is left for future investigations. 
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